From: wolf Date: Fri, 27 Apr 2001 22:01:17 +0000 (+0000) Subject: Add filtered matrix, test case, doc updates, etc. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b9572da587ebd15d1d32d6c153c5c7a01cbd48b6;p=dealii-svn.git Add filtered matrix, test case, doc updates, etc. git-svn-id: https://svn.dealii.org/trunk@4498 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index a7a669649c..ef34ba6759 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -502,24 +502,27 @@ class MatrixCreator * * @sect3{Boundary conditions} * - * The @p{apply_boundary_values} function inserts boundary conditions of - * into a system of equations. To actually do this you have to specify - * a list of degree of freedom indices along with the values these degrees of - * freedom shall assume. To see how to get such a list, see the discussion - * of the @ref{VectorTools}@p{::interpolate_boundary_values} function. - * - * The inclusion into the assemblage process is as follows: when the matrix and - * vectors are set up, a list of nodes subject to dirichlet bc is made and - * matrix and vectors are changed accordingly. This is done by deleting all - * entries in the matrix in the line of this degree of freedom, setting the - * main diagonal entry to one and the right hand side element to the - * boundary value at this node. This forces this node's value to be as specified. - * To decouple the remaining linear system of equations and to make the system + * The @p{apply_boundary_values} function inserts boundary conditions + * of into a system of equations. To actually do this you have to + * specify a list of degree of freedom indices along with the values + * these degrees of freedom shall assume. To see how to get such a + * list, see the discussion of the + * @ref{VectorTools}@p{::interpolate_boundary_values} function. + * + * The inclusion into the assemblage process is as follows: when the + * matrix and vectors are set up, a list of nodes subject to dirichlet + * bc is made and matrix and vectors are changed accordingly. This is + * done by deleting all entries in the matrix in the line of this + * degree of freedom, setting the main diagonal entry to one and the + * right hand side element to the boundary value at this node. This + * forces this node's value to be as specified. To decouple the + * remaining linear system of equations and to make the system * symmetric again (at least if it was before), one Gauss elimination - * step is performed with this line, by adding this (now almost empty) line to - * all other lines which couple with the given degree of freedom and thus - * eliminating all coupling between this degree of freedom and others. Now - * also the column consists only of zeroes, apart from the main diagonal entry. + * step is performed with this line, by adding this (now almost empty) + * line to all other lines which couple with the given degree of + * freedom and thus eliminating all coupling between this degree of + * freedom and others. Now also the column consists only of zeroes, + * apart from the main diagonal entry. * * Finding which rows contain an entry in the column for which we are * presently performing a Gauss elimination step is either difficult @@ -558,14 +561,15 @@ class MatrixCreator * would be @p{O(N*sqrt(N)*log(m))} for the general case; the latter * is too expensive to be performed. * - * It seems as if we had to make clear not to overwrite the lines of other - * boundary nodes when doing the Gauss elimination step. However, since we - * reset the right hand side when passing such a node, it is not a problem - * to change the right hand side values of other boundary nodes not yet - * processed. It would be a problem to change those entries of nodes already - * processed, but since the matrix entry of the present column on the row - * of an already processed node is zero, the Gauss step does not change - * the right hand side. We need therefore not take special care of other + * It seems as if we had to make clear not to overwrite the lines of + * other boundary nodes when doing the Gauss elimination + * step. However, since we reset the right hand side when passing such + * a node, it is not a problem to change the right hand side values of + * other boundary nodes not yet processed. It would be a problem to + * change those entries of nodes already processed, but since the + * matrix entry of the present column on the row of an already + * processed node is zero, the Gauss step does not change the right + * hand side. We need therefore not take special care of other * boundary nodes. * * To make solving faster, we preset the solution vector with the @@ -585,6 +589,17 @@ class MatrixCreator * set the entry to the mean of the other diagonal entries, but this * seems to be too expensive. * + * In some cases, it might be interesting to solve several times with + * the same matrix, but for different right hand sides or boundary + * values. However, since the modification for boundary values of the + * right hand side vector depends on the original matrix, this is not + * possible without storing the original matrix somewhere and applying + * the @p{apply_boundary_conditions} function to a copy of it each + * time we want to solve. In that case, you can use the + * @ref{FilteredMatrix} class in the @p{LAC} sublibrary. There you can + * also find a formal (mathematical) description of the process of + * modifying the matrix and right hand side vectors for boundary + * values. * * @author Wolfgang Bangerth, 1998, 2000 */ @@ -597,6 +612,11 @@ class MatrixTools : public MatrixCreator * to the system matrix and vectors * as described in the general * documentation. + * + * For a replacement function, + * see the documentation of the + * @ref{FilteredMatrix} class in + * the @p{LAC} sublibrary. */ template static void @@ -614,6 +634,11 @@ class MatrixTools : public MatrixCreator * documentation. This function * works for block sparse * matrices and block vectors + * + * For a replacement function, + * see the documentation of the + * @ref{FilteredMatrix} class in + * the @p{LAC} sublibrary. */ static void apply_boundary_values (const std::map &boundary_values, diff --git a/deal.II/doc/future.html b/deal.II/doc/future.html index cf74def2ff..807c0a9151 100644 --- a/deal.II/doc/future.html +++ b/deal.II/doc/future.html @@ -130,8 +130,9 @@ and SolutionTransfer classes, the AssertThrow macro, 1d programs, curved boundaries and different mappings, using the - MatrixTools::create_*_matrix functions, - and output of more than one variable). This is certainly + MatrixTools::create_*_matrix + functions, using the FilteredMatrix + class, and output of more than one variable). This is certainly something that should be improved, but is rather time consuming.

diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 30f2e8e513..6e104a6158 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -249,6 +249,17 @@ documentation, etc. (WB 2001/04/25)

+
  • + New: The FilteredMatrix class is a + replacement for the MatrixTools::apply_boundary_values + function for cases where you would like to solve several times + with the same matrix, either for different right hand sides, or + for different boundary values. +
    + (WB 2001/04/27) +

    +
  • New: There is now a function Vector::scale(Vector) diff --git a/deal.II/lac/include/lac/filtered_matrix.h b/deal.II/lac/include/lac/filtered_matrix.h new file mode 100644 index 0000000000..03dc7f33a4 --- /dev/null +++ b/deal.II/lac/include/lac/filtered_matrix.h @@ -0,0 +1,629 @@ +//---------------------------- filtered_matrix.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2001 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. +// +//---------------------------- filtered_matrix.h --------------------------- +#ifndef __deal2__filtered_matrix_h +#define __deal2__filtered_matrix_h + + + +#include +#include +#include +#include + + +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 + * 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 + * degrees of freedom (e.g. for Dirichlet boundary nodes, thus the + * index $D$) and zeroes for all other diagonal entries, and $g_D$ + * having the requested nodal values for these constrained + * nodes. Thus, the underdetermined equation $B_D x = g_D$ fixes only + * the constrained nodes and does not impose any condition on the + * others. We note that $B_D B_D^T = 1_D$, where $1_D$ is the identity + * matrix with dimension as large as the number of constrained degrees + * of freedom. Likewise, $B_D^T B_D$ is the diagonal matrix with + * diagonal entries $0$ or $1$ that, when applied to a vector, leaves + * all constrained nodes untouched and deletes all unconstrained ones. + * + * For solving such a system of equations, we first write down the + * Lagrangian $L=1/2 x^T A x - x^T b + l^T B_D x$, where $l$ + * is a Lagrange multiplier for the constraints. The stationarity + * condition then reads + * \begin{verbatim} + * [ A B_D^T ] [x] = [b ] + * [ B_D 0 ] [l] = [g_D] + * \begin{verbatim} + * + * The first equation then reads $B_D^T l = b-Ax$. On the other hand, + * if we left-multiply the first equation by $B_D^T B_D$, we obtain + * $B_D^T B_D A x + B_D^T l = B_D^T B_D b$ after equating $B_D B_D^T$ + * to the identity matrix. Inserting the previous equality, this + * yields $(A - B_D^T B_D A) x = (1 - B_D^T B_D)b$. Since + * $x=(1 - B_D^T B_D) x + B_D^T B_D x = (1 - B_D^T B_D) x + B_D^T g_D$, + * we can restate the linear system: + * $A_D x = (1 - B_D^T B_D)b - (1 - B_D^T B_D) A B^T g_D$, where + * $A_D = (1 - B_D^T B_D) A (1 - B_D^T B_D)$ is the matrix where all + * rows and columns corresponding to constrained nodes have been deleted. + * + * The last system of equation only defines the value of the + * unconstrained nodes, while the constrained ones are determined by + * the equation $B_D x = g_D$. We can combine these two linear systems + * by using the zeroed out rows of $A_D$: if we set the diagonal to + * $1$ and the corresponding zeroed out element of the right hand side + * to that of $g_D$, then this fixes the constrained elements as + * well. We can write this as follows: + * $A_X x = (1 - B_D^T B_D)b - (1 - B_D^T B_D) A B^T g_D + B_D^T g_D$, + * where $A_X = A_D + B_D^T B_D$. Note that the two parts of the + * latter matrix operate on disjoint subspaces (the first on the + * unconstrained nodes, the latter on the constrained ones). + * + * In iterative solvers, it is not actually necessary to compute $A_X$ + * explicitely, since only matrix-vector operations need to be + * performed. This can be done in a three-step procedure that first + * clears all elements in the incoming vector that belong to + * constrained nodes, then performs the product with the matrix $A$, + * then clears again. This class is a wrapper to this procedure, it + * takes a pointer to a matrix with which to perform matrix-vector + * products, and does the cleaning of constrained elements itself. + * This class therefore implements an overloaded @p{vmult} function + * that does the matrix-vector product, as well as @p{Tvmult} for + * transpose matrix-vector multiplication and @p{residual} for + * residual computation, and can thus be used as a matrix replacement + * in lineaer solvers. + * + * It also has the ability to generate the modification of the right + * hand side, through the @ref{apply_constraints} function. + * + * + * @sect3{Connection to other classes} + * + * The function @p{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 + * @p{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. + * + * + * @sect3{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 + * @ref{set_referenced_matrix} function), specify the list of boundary + * values or other constraints (through the @ref{add_constraints} + * function), and then for each required solution modify the right + * hand side vector (through @ref{apply_constraints}) and use this + * object as matrix object in a linear solver. As linear solvers + * should only use @ref{vmult} and @ref{residual} functions of a + * matrix class, this class should be as a good a matrix as any other + * for that purpose. + * + * Furthermore, also the @ref{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: + * @begin{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); + * @end{verbatim} + * + * + * @sect3{Template arguments} + * + * This class takes as template arguments a matrix and a vector + * class. The former must provide @p{vmult}, @p{Tvmult}, and + * @p{residual} member function that operate on the vector type (the + * second template argument). The latter template parameter must + * provide access to indivual elements through @p{operator()}, + * assignment through @p{operator=}. + * + * + * @sect3{Thread-safety} + * + * The functions that operate as a matrix and do not change the + * internal state of this object are synchronised and thus + * threadsafe. You need not serialize calls to @p{vmult} or + * @p{residual} therefore. Because these functions require the use of + * a temporary, they block mutual execution, however. It is necessary + * to allocate this temporary vector in class space since otherwise we + * would have to allocate such a vector each time one of the member + * functions is called (which may be very often for slowly converging + * linear systems), which would be a serious performance + * bottleneck. If you don't want this serialization of operations, you + * have to use several objects of this type. + * + * @author Wolfgang Bangerth 2001 + */ +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 typename std::pair IndexValuePair; + + /** + * Default constructor. You will + * have to set the matrix to be + * used later using the + * @{set_referenced_matrix} + * function. + */ + FilteredMatrix (); + + /** + * Copy constructor. Use the + * matrix and the constraints set + * in the given object for the + * present one as well. + */ + FilteredMatrix (const FilteredMatrix &fm); + + /** + * Constructor. Use the given + * matrix for future operations. + */ + FilteredMatrix (const Matrix &matrix); + + /** + * Copy operator. Take over + * matrix and constraints from + * the other object. + */ + FilteredMatrix & operator = (const FilteredMatrix &fm); + + /** + * Set the matrix to be used + * further on. You will probably + * also want to call the + * @ref{clear_constraints} + * function if constraits were + * previously added. + */ + void set_referenced_matrix (const Matrix &m); + + /** + * 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 + * this object. The actual data + * type of this list must be so + * that dereferenced iterators + * are pairs of indices and the + * corresponding values to be + * enforced on the respective + * solution vector's entry. Thus, + * the data type might be, for + * example, a @{std::list} or + * @p{std::vector} of + * @ref{IndexValuePair} objects, + * but also a + * @p{std::map}. + * + * It is an error if the argument + * contains an entry for a degree + * of freedom that has already + * been constrained + * previously. Furthermore, it is + * assumed that the list of + * constraints is sorted. If the + * input argument is a + * @p{std::map}, + * this is automatically the + * case. + */ + template + void add_constraints (const ConstraintList &new_constraints); + + /** + * Delete the list of constraints + * presently in use. + */ + void clear_constraints (); + + /** + * Apply the constraints to a + * right hand side vector. This + * needs to be done before + * starting to solve with the + * filtered matrix. If the matrix + * is symmetric, set the second + * parameter to @p{true} to use a + * faster algorithm. + */ + 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$ + * being this matrix. (This + * matrix is the filtered one to + * which we store a reference.) + */ + void vmult (Vector &dst, + const Vector &src) const; + + /** + * Matrix-vector multiplication: + * let $dst = M^T*src$ with $M$ + * being this matrix. This + * function does the same as + * @p{vmult} but takes the + * transposed matrix. (This + * matrix is the filtered one to + * which we store a reference.) + * + * Because we need to use a + * temporary variable and since + * we only allocate that each + * time the matrix changed, this + * function only works for square + * matrices. + */ + void Tvmult (Vector &dst, + 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 @p{Mx=b}, where the + * residual is defined to be + * @p{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. + * + * Note that it is assumed that + * @{b} is a vector that has been + * treated by the + * @ref{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. + */ + value_type residual (Vector &dst, + const Vector &x, + const Vector &b) 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}. + */ + void precondition_Jacobi (Vector &dst, + const Vector &src, + const value_type omega = 1.) const; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. Since we are + * not the owner of the matrix + * referenced, its memory + * consumption is not included. + */ + unsigned int memory_consumption () const; + + private: + /** + * Declare an abbreviation for an + * iterator into the array + * constraint pairs, since that + * data type is so often used and + * is rather awkward to write out + * each time. + */ + typedef typename std::vector::const_iterator const_index_value_iterator; + + /** + * Helper class used to sort + * pairs of indices and + * values. Only the index is + * considered as sort key. + */ + struct PairComparison + { + /** + * Function comparing the + * pairs @p{i1} and @p{i2} + * for their keys. + */ + bool operator () (const IndexValuePair &i1, + const IndexValuePair &i2) const; + }; + + /** + * Pointer to the sparsity + * pattern used for this + * matrix. In order to guarantee + * that it is not deleted while + * still in use, we subscribe to + * it using the @p{SmartPointer} + * class. + */ + SmartPointer matrix; + + /** + * Sorted list of pairs denoting + * the index of the variable and + * the value to which it shall be + * fixed. + */ + std::vector constraints; + + /** + * Vector to be used as temporary + * storage. Since memory + * allocation is expensive, we do + * not want to allocate temporary + * vectors in each call to + * matrix-vector function, so we + * rather allocate it only once + * and then reuse it over and + * over again. Note that in a + * multithreaded environment, we + * have to synchronise access to + * this vector. + */ + mutable Vector tmp_vector; + + /** + * Mutex used to synchronise use + * of the temporary vector. + */ + mutable Threads::ThreadMutex tmp_mutex; + + /** + * Do the pre-filtering step, + * i.e. zero out those components + * that belong to constrained + * degrees of freedom. + */ + void pre_filter (Vector &v) const; + + /** + * Do the postfiltering step, + * i.e. set constrained degrees + * of freedom to the value of the + * input vector, as the matrix + * contains only ones on the + * diagonal for these degrees of + * freedom. + */ + 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. + */ + void allocate_tmp_vector (); + + /** + * Determine all entries in the + * given column of the matrix + * except for the diagonal entry + * and return their index/value + * pairs. If the matrix is + * symmetric, use a faster + * algorithm. + * + * This function needs to be + * specialised for the different + * matrix types. + */ + void get_column_entries (const unsigned int index, + std::vector &column_entries, + const bool matrix_is_symmetric) const; +}; + + +/*---------------------- Inline functions -----------------------------------*/ + + +template +inline +bool +FilteredMatrix::PairComparison:: +operator () (const IndexValuePair &i1, + const IndexValuePair &i2) const +{ + return (i1.first < i2.first); +}; + + + +template +template +void +FilteredMatrix:: +add_constraints (const ConstraintList &new_constraints) +{ + // add new constraints to end + const unsigned int old_size = constraints.size(); + constraints.reserve (old_size + new_constraints.size()); + constraints.insert (constraints.end(), + new_constraints.begin(), + new_constraints.end()); + // then merge the two arrays to + // form one sorted one + std::inplace_merge (constraints.begin(), + constraints.begin()+old_size, + constraints.end(), + PairComparison()); +//TODO:[WB] Use equal_range etc to assert that the array is indeed sorted +}; + + + +template +inline +const Matrix & +FilteredMatrix::get_referenced_matrix () const +{ + return *matrix; +}; + + + +template +inline +unsigned int FilteredMatrix::m () const +{ + return matrix->m(); +}; + + + +template +inline +unsigned int FilteredMatrix::n () const +{ + return matrix->n(); +}; + + + +/*---------------------------- filtered_matrix.h ---------------------------*/ + +#endif +/*---------------------------- filtered_matrix.h ---------------------------*/ + + diff --git a/deal.II/lac/include/lac/filtered_matrix.templates.h b/deal.II/lac/include/lac/filtered_matrix.templates.h new file mode 100644 index 0000000000..f24a79bda1 --- /dev/null +++ b/deal.II/lac/include/lac/filtered_matrix.templates.h @@ -0,0 +1,424 @@ +//---------------------------- filtered_matrix.templates.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001 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. +// +//---------------------------- filtered_matrix.templates.h --------------------------- +#ifndef __deal2__filtered_matrix_templates_h +#define __deal2__filtered_matrix_templates_h + + +#include +#include +#include +#include +#include +#include + + +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 +{ + // array that will hold the pairs + // of index/value of all nonzero + // entries in a given column + std::vector column_entries; + + // iterate over all constraints and + // treat them one after the other + const_index_value_iterator i = constraints.begin(); + const const_index_value_iterator e = constraints.end(); + for (; i!=e; ++i) + { + // define abbreviations + const unsigned index = i->first; + const value_type value = i->second; + + // check whether the value is + // zero, since in that case we do + // not have to modify other nodes + if (value != 0) + { + // first clear array of + // previous content + column_entries.clear (); + + // then get all entries in + // the present column + get_column_entries (index, column_entries, matrix_is_symmetric); + + // modify rhs for each entry + const_index_value_iterator col = column_entries.begin(); + const const_index_value_iterator col_end = column_entries.end(); + for (; col!=col_end; ++col) + v(col->first) -= col->second * value; + }; + }; + + + // finally set constrained + // entries themselves. we can't + // do it in the above loop + // since we might end up + // modifying an entry that we + // have already set if + // constrained dofs couple to + // each other + for (i=constraints.begin(); i!=e; ++i) + v(i->first) = i->second; +}; + + + +template <> +void +FilteredMatrix,Vector >:: +get_column_entries (const unsigned int index, + std::vector &column_entries, + const bool matrix_is_symmetric) const +{ + // depending on whether the matrix + // can be assumed symmetric or not, + // either use a fast or a slow + // algorithm + if (matrix_is_symmetric == true) + // ok, matrix is symmetric. we + // may determine the matrix + // entries in this column by + // looking at the matrix entries + // in this row which is + // significantly faster since we + // can traverse them linearly and + // do not have to check each row + // for the possible existence of + // a matrix entry + { + const unsigned int * + col_nums = &(matrix->get_sparsity_pattern().get_column_numbers() + [matrix->get_sparsity_pattern().get_rowstart_indices()[index]]); + const unsigned int + row_length = matrix->get_sparsity_pattern().row_length(index); + + for (unsigned int i=0; iget_sparsity_pattern()(row,index); + if (global_index != SparsityPattern::invalid_entry) + column_entries.push_back (std::make_pair(row, + (*matrix)(row,index))); + }; + }; +}; + + + +template <> +void +FilteredMatrix,BlockVector >:: +get_column_entries (const unsigned int /*index*/, + std::vector &/*column_entries*/, + const bool /*matrix_is_symmetric*/) const +{ + // presently not implemented, but + // should be fairly simple to do + Assert (false, ExcNotImplemented()); +}; + + + +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)); +}; + + + +template <> +void +FilteredMatrix,Vector >:: +allocate_tmp_vector () +{ + tmp_mutex.acquire (); + tmp_vector.reinit (matrix->n()); + tmp_mutex.release (); +}; + + + +template <> +void +FilteredMatrix,Vector >:: +allocate_tmp_vector () +{ + tmp_mutex.acquire (); + tmp_vector.reinit (matrix->n()); + tmp_mutex.release (); +}; + + + +template <> +void +FilteredMatrix,BlockVector >:: +allocate_tmp_vector () +{ + tmp_mutex.acquire (); + tmp_vector.reinit (matrix->n()); + tmp_mutex.release (); +}; + + + +template <> +void +FilteredMatrix,BlockVector >:: +allocate_tmp_vector () +{ + tmp_mutex.acquire (); + tmp_vector.reinit (matrix->n()); + tmp_mutex.release (); +}; + + +#endif diff --git a/deal.II/lac/source/filtered_matrix.cc b/deal.II/lac/source/filtered_matrix.cc new file mode 100644 index 0000000000..40dc8c4242 --- /dev/null +++ b/deal.II/lac/source/filtered_matrix.cc @@ -0,0 +1,19 @@ +//---------------------------- filtered_matrix.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2001 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. +// +//---------------------------- filtered_matrix.cc --------------------------- + + +#include + + +template class FilteredMatrix,Vector >; +template class FilteredMatrix,BlockVector >; diff --git a/tests/deal.II/Makefile.in b/tests/deal.II/Makefile.in index 2f76a715b2..cbd9661400 100644 --- a/tests/deal.II/Makefile.in +++ b/tests/deal.II/Makefile.in @@ -40,13 +40,13 @@ mglocal.exe : mglocal.go $(lib-2d) second_derivatives.exe : second_derivatives.go $(lib-2d) $(libraries) wave-test-3.exe : wave-test-3.go $(lib-2d) $(libraries) support_point_map.exe : support_point_map.go $(lib-1d) $(lib-2d) $(lib-3d) $(libraries) - +filtered_matrix.exe : filtered_matrix.go $(lib-1d) $(lib-2d) $(lib-3d) $(libraries) tests = grid_test dof_test data_out derivatives gradients constraints mg \ mglocal block_matrices second_derivatives derivative_approximation \ matrices error_estimator intergrid_constraints intergrid_map \ - wave-test-3 dof_renumbering support_point_map + wave-test-3 dof_renumbering support_point_map filtered_matrix ############################################################ diff --git a/tests/deal.II/filtered_matrix.cc b/tests/deal.II/filtered_matrix.cc new file mode 100644 index 0000000000..7864d44025 --- /dev/null +++ b/tests/deal.II/filtered_matrix.cc @@ -0,0 +1,214 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// (c) 2001 the deal.II authors +// +// purpose of this test: +// +// compare results with boundary values eliminated from matrix and +// vector, and with boundary values treated by filtering +// +// +// Method: +// +// solve (u,v) = (f,v) +// +//---------------------------------------------------------------------- + + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + + +void +solve_filtered (std::map &bv, + SparseMatrix &A, + Vector &u, + Vector &f) +{ + FilteredMatrix > A1 (A); + A1.add_constraints (bv); + + SolverControl control (1000, 1.e-10, false, false); + PrimitiveVectorMemory > mem; + SolverCG > solver (control, mem); + PreconditionJacobi > > prec; + prec.initialize (A1, 1.2); + + Vector f1 (f.size()); + f1 = f; + A1.apply_constraints (f1, true); + + solver.solve (A1, u, f1, prec); + + for (typename std::map::const_iterator i=bv.begin(); + i!=bv.end(); ++i) + Assert (std::fabs(u(i->first) - i->second) < 1e-8, + ExcInternalError()); +}; + + + +template +void +solve_eliminated (std::map &bv, + SparseMatrix &A, + Vector &u, + Vector &f) +{ + MatrixTools::apply_boundary_values (bv, A, u, f); + + SolverControl control (1000, 1.e-10, false, false); + PrimitiveVectorMemory > mem; + SolverCG > solver (control, mem); + PreconditionJacobi<> prec; + prec.initialize (A, 1.2); + + solver.solve (A, u, f, prec); +}; + + + +template +void +check () +{ + Triangulation tr; + + CosineFunction cosine; + + if (dim==2) + GridGenerator::hyper_ball(tr); + else + GridGenerator::hyper_cube(tr, -1,1); + + tr.refine_global (5-dim); + + MappingQ mapping(2); + FE_Q element(1); + QGauss4 quadrature; + + DoFHandler dof(tr); + dof.distribute_dofs(element); + + FEValues fe (mapping, element, quadrature, + update_values | update_gradients + | update_q_points | update_JxW_values); + + vector global_dofs (element.dofs_per_cell); + vector function (quadrature.n_quadrature_points); + + Vector f (dof.n_dofs ()); + + SparsityPattern A_pattern (dof.n_dofs (), dof.n_dofs (), + dof.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern(dof, A_pattern); + A_pattern.compress (); + + SparseMatrix A(A_pattern); + + DoFHandler::active_cell_iterator cell = dof.begin_active(); + const DoFHandler::cell_iterator end = dof.end(); + + for (; cell != end;++cell) + { + fe.reinit(cell); + cell->get_dof_indices (global_dofs); + cosine.value_list (fe.get_quadrature_points(), function); + + for (unsigned int k=0;k grad_v = fe.shape_grad(i,k); + + double rhs = dx * v * (function[k]); + + f(global_dofs[i]) += rhs; + for (unsigned int j=0;j grad_u = fe.shape_grad (j,k); + double el = dx * (grad_u*grad_v); + A.add (global_dofs[i], global_dofs[j], el); + } + } + } + } + + // interpolate boundary values + std::map bv; + VectorTools::interpolate_boundary_values (mapping, dof, 0, cosine, bv); + // the cosine has too many zero + // values on the boundary of the + // domain, so reset the elements to + // some other value + for (typename std::map::iterator i=bv.begin(); + i!=bv.end(); ++i) + i->second = std::sin(i->second+0.5)+1.0; + + // first solve filtered. this does + // not change the matrix + Vector u_filtered (dof.n_dofs ()); + solve_filtered (bv, A, u_filtered, f); + + // then solve by eliminating in the + // matrix. since this changes the + // matrix, this call must come + // second + Vector u_eliminated (dof.n_dofs ()); + solve_eliminated (bv, A, u_eliminated, f); + + // output and check + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/filtered_matrix.checked b/tests/deal.II/filtered_matrix.checked new file mode 100644 index 0000000000..68f04927b2 --- /dev/null +++ b/tests/deal.II/filtered_matrix.checked @@ -0,0 +1,480 @@ + +DEAL:1d::1.48 +DEAL:1d::1.64 +DEAL:1d::1.79 +DEAL:1d::1.94 +DEAL:1d::2.08 +DEAL:1d::2.21 +DEAL:1d::2.33 +DEAL:1d::2.43 +DEAL:1d::2.52 +DEAL:1d::2.59 +DEAL:1d::2.65 +DEAL:1d::2.69 +DEAL:1d::2.72 +DEAL:1d::2.74 +DEAL:1d::2.75 +DEAL:1d::2.75 +DEAL:1d::2.75 +DEAL:2d::1.64 +DEAL:2d::1.72 +DEAL:2d::1.74 +DEAL:2d::1.68 +DEAL:2d::1.77 +DEAL:2d::1.78 +DEAL:2d::1.80 +DEAL:2d::1.77 +DEAL:2d::1.72 +DEAL:2d::1.80 +DEAL:2d::1.81 +DEAL:2d::1.81 +DEAL:2d::1.82 +DEAL:2d::1.83 +DEAL:2d::1.82 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.83 +DEAL:2d::1.79 +DEAL:2d::1.75 +DEAL:2d::1.81 +DEAL:2d::1.78 +DEAL:2d::1.80 +DEAL:2d::1.81 +DEAL:2d::1.77 +DEAL:2d::1.78 +DEAL:2d::1.80 +DEAL:2d::1.82 +DEAL:2d::1.72 +DEAL:2d::1.74 +DEAL:2d::1.64 +DEAL:2d::1.68 +DEAL:2d::1.72 +DEAL:2d::1.77 +DEAL:2d::1.79 +DEAL:2d::1.81 +DEAL:2d::1.75 +DEAL:2d::1.78 +DEAL:2d::1.81 +DEAL:2d::1.83 +DEAL:2d::1.83 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.86 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.83 +DEAL:2d::1.80 +DEAL:2d::1.82 +DEAL:2d::1.84 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.84 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.88 +DEAL:2d::1.88 +DEAL:2d::1.83 +DEAL:2d::1.80 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.82 +DEAL:2d::1.85 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.88 +DEAL:2d::1.87 +DEAL:2d::1.86 +DEAL:2d::1.84 +DEAL:2d::1.87 +DEAL:2d::1.86 +DEAL:2d::1.74 +DEAL:2d::1.72 +DEAL:2d::1.77 +DEAL:2d::1.80 +DEAL:2d::1.78 +DEAL:2d::1.77 +DEAL:2d::1.79 +DEAL:2d::1.81 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.83 +DEAL:2d::1.82 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.80 +DEAL:2d::1.82 +DEAL:2d::1.81 +DEAL:2d::1.83 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.86 +DEAL:2d::1.88 +DEAL:2d::1.88 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.85 +DEAL:2d::1.86 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.86 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.83 +DEAL:2d::1.87 +DEAL:2d::1.88 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.86 +DEAL:2d::1.84 +DEAL:2d::1.87 +DEAL:2d::1.86 +DEAL:2d::1.84 +DEAL:2d::1.82 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.80 +DEAL:2d::1.78 +DEAL:2d::1.81 +DEAL:2d::1.80 +DEAL:2d::1.82 +DEAL:2d::1.80 +DEAL:2d::1.78 +DEAL:2d::1.77 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.79 +DEAL:2d::1.77 +DEAL:2d::1.75 +DEAL:2d::1.72 +DEAL:2d::1.74 +DEAL:2d::1.72 +DEAL:2d::1.68 +DEAL:2d::1.64 +DEAL:2d::1.88 +DEAL:2d::1.88 +DEAL:2d::1.89 +DEAL:2d::1.88 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.90 +DEAL:2d::1.89 +DEAL:2d::1.90 +DEAL:2d::1.89 +DEAL:2d::1.90 +DEAL:2d::1.90 +DEAL:2d::1.90 +DEAL:2d::1.90 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.88 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.88 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.88 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.88 +DEAL:2d::1.88 +DEAL:2d::1.89 +DEAL:2d::1.90 +DEAL:2d::1.90 +DEAL:2d::1.90 +DEAL:2d::1.90 +DEAL:2d::1.90 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.90 +DEAL:2d::1.89 +DEAL:2d::1.88 +DEAL:2d::1.87 +DEAL:2d::1.88 +DEAL:2d::1.88 +DEAL:2d::1.88 +DEAL:2d::1.87 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.88 +DEAL:2d::1.88 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.88 +DEAL:2d::1.90 +DEAL:2d::1.89 +DEAL:2d::1.89 +DEAL:2d::1.88 +DEAL:2d::1.88 +DEAL:2d::1.87 +DEAL:2d::1.88 +DEAL:2d::1.87 +DEAL:2d::1.72 +DEAL:2d::1.74 +DEAL:2d::1.77 +DEAL:2d::1.78 +DEAL:2d::1.80 +DEAL:2d::1.77 +DEAL:2d::1.80 +DEAL:2d::1.81 +DEAL:2d::1.81 +DEAL:2d::1.82 +DEAL:2d::1.83 +DEAL:2d::1.82 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.83 +DEAL:2d::1.79 +DEAL:2d::1.81 +DEAL:2d::1.80 +DEAL:2d::1.81 +DEAL:2d::1.77 +DEAL:2d::1.78 +DEAL:2d::1.80 +DEAL:2d::1.82 +DEAL:2d::1.72 +DEAL:2d::1.74 +DEAL:2d::1.64 +DEAL:2d::1.68 +DEAL:2d::1.72 +DEAL:2d::1.77 +DEAL:2d::1.79 +DEAL:2d::1.81 +DEAL:2d::1.75 +DEAL:2d::1.78 +DEAL:2d::1.81 +DEAL:2d::1.83 +DEAL:2d::1.83 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.86 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.83 +DEAL:2d::1.80 +DEAL:2d::1.82 +DEAL:2d::1.84 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.84 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.83 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.86 +DEAL:2d::1.74 +DEAL:2d::1.72 +DEAL:2d::1.77 +DEAL:2d::1.80 +DEAL:2d::1.78 +DEAL:2d::1.77 +DEAL:2d::1.79 +DEAL:2d::1.81 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.83 +DEAL:2d::1.82 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.80 +DEAL:2d::1.82 +DEAL:2d::1.81 +DEAL:2d::1.83 +DEAL:2d::1.84 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.86 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.85 +DEAL:2d::1.86 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.86 +DEAL:2d::1.85 +DEAL:2d::1.84 +DEAL:2d::1.83 +DEAL:2d::1.87 +DEAL:2d::1.87 +DEAL:2d::1.86 +DEAL:2d::1.84 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.81 +DEAL:2d::1.80 +DEAL:2d::1.82 +DEAL:2d::1.80 +DEAL:2d::1.78 +DEAL:2d::1.77 +DEAL:2d::1.83 +DEAL:2d::1.81 +DEAL:2d::1.79 +DEAL:2d::1.77 +DEAL:2d::1.74 +DEAL:2d::1.72 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.53 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.59 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.59 +DEAL:3d::1.63 +DEAL:3d::1.59 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.53 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.59 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.53 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.55 +DEAL:3d::1.59 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.53 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.53 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.55 +DEAL:3d::1.59 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.53 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.53 +DEAL:3d::1.55 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.53 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48 +DEAL:3d::1.48