*
* @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
* 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
* 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
*/
* 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 <typename number>
static void
* 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<unsigned int,double> &boundary_values,
and <code class="class">SolutionTransfer</code> classes, the
<code class="class">AssertThrow</code> macro, 1d programs,
curved boundaries and different mappings, using the
- <code class="member">MatrixTools::create_*_matrix</code> functions,
- and output of more than one variable). This is certainly
+ <code class="member">MatrixTools::create_*_matrix</code>
+ functions, using the <code class="class">FilteredMatrix</code>
+ class, and output of more than one variable). This is certainly
something that should be improved, but is rather time
consuming.
</p>
(WB 2001/04/25)
</p>
+ <li> <p>
+ New: The <code class="class">FilteredMatrix</code> class is a
+ replacement for the <code
+ class="class">MatrixTools::apply_boundary_values</code>
+ 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.
+ <br>
+ (WB 2001/04/27)
+ </p>
+
<li> <p>
New: There is now a function <code
class="member">Vector::scale(Vector)</code>
--- /dev/null
+//---------------------------- 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 <base/smartpointer.h>
+#include <base/thread_management.h>
+#include <vector>
+#include <algorithm>
+
+
+template <typename number> 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<SparseMatrix<double> > filtered_A (A);
+ * filtered_A.add_constraints (boundary_values);
+ *
+ * // set up a linear solver
+ * SolverControl control (1000, 1.e-10, false, false);
+ * PrimitiveVectorMemory<Vector<double> > mem;
+ * SolverCG<Vector<double> > solver (control, mem);
+ *
+ * // set up a preconditioner object
+ * PreconditionJacobi<FilteredMatrix<SparseMatrix<double> > > 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 Matrix, class Vector=::Vector<typename Matrix::value_type> >
+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<unsigned int,value_type> 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<unsigned,value_type>}.
+ *
+ * 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<unsigned,value_type>},
+ * this is automatically the
+ * case.
+ */
+ template <class ConstraintList>
+ 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<IndexValuePair>::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<const Matrix> matrix;
+
+ /**
+ * Sorted list of pairs denoting
+ * the index of the variable and
+ * the value to which it shall be
+ * fixed.
+ */
+ std::vector<IndexValuePair> 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<IndexValuePair> &column_entries,
+ const bool matrix_is_symmetric) const;
+};
+
+
+/*---------------------- Inline functions -----------------------------------*/
+
+
+template <class Matrix, class Vector>
+inline
+bool
+FilteredMatrix<Matrix,Vector>::PairComparison::
+operator () (const IndexValuePair &i1,
+ const IndexValuePair &i2) const
+{
+ return (i1.first < i2.first);
+};
+
+
+
+template <class Matrix, class Vector>
+template <class ConstraintList>
+void
+FilteredMatrix<Matrix,Vector>::
+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 <class Matrix, class Vector>
+inline
+const Matrix &
+FilteredMatrix<Matrix,Vector>::get_referenced_matrix () const
+{
+ return *matrix;
+};
+
+
+
+template <class Matrix, class Vector>
+inline
+unsigned int FilteredMatrix<Matrix,Vector>::m () const
+{
+ return matrix->m();
+};
+
+
+
+template <class Matrix, class Vector>
+inline
+unsigned int FilteredMatrix<Matrix,Vector>::n () const
+{
+ return matrix->n();
+};
+
+
+
+/*---------------------------- filtered_matrix.h ---------------------------*/
+
+#endif
+/*---------------------------- filtered_matrix.h ---------------------------*/
+
+
--- /dev/null
+//---------------------------- 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 <base/memory_consumption.h>
+#include <lac/filtered_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+
+
+template <class Matrix, class Vector>
+FilteredMatrix<Matrix,Vector>::
+FilteredMatrix ()
+{};
+
+
+
+template <class Matrix, class Vector>
+FilteredMatrix<Matrix,Vector>::
+FilteredMatrix (const FilteredMatrix &fm)
+ :
+ Subscriptor (),
+ constraints (fm.constraints)
+{
+ set_referenced_matrix (*fm.matrix);
+};
+
+
+
+template <class Matrix, class Vector>
+FilteredMatrix<Matrix,Vector>::
+FilteredMatrix (const Matrix &m)
+{
+ set_referenced_matrix (m);
+};
+
+
+
+template <class Matrix, class Vector>
+FilteredMatrix<Matrix,Vector> &
+FilteredMatrix<Matrix,Vector>::operator = (const FilteredMatrix &fm)
+{
+ set_referenced_matrix (*fm.matrix);
+ constraints = fm.constraints;
+ return *this;
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::
+set_referenced_matrix (const Matrix &m)
+{
+ matrix = &m;
+ allocate_tmp_vector ();
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::clear_constraints ()
+{
+ // swap vectors to release memory
+ std::vector<IndexValuePair> empty;
+ constraints.swap (empty);
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::
+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<IndexValuePair> 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<SparseMatrix<double>,Vector<double> >::
+get_column_entries (const unsigned int index,
+ std::vector<IndexValuePair> &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; i<row_length; ++i)
+ {
+ const unsigned int c = *(col_nums+i);
+
+ // if not diagonal entry,
+ // add to list
+ if (c != index)
+ column_entries.push_back (std::make_pair(c, (*matrix)(c,index)));
+ };
+ }
+ else
+ {
+ // otherwise check each row for
+ // occurrence of an entry in
+ // this column
+ for (unsigned int row=0; row<n(); ++row)
+ if (row != index)
+ {
+ const unsigned int
+ global_index = matrix->get_sparsity_pattern()(row,index);
+ if (global_index != SparsityPattern::invalid_entry)
+ column_entries.push_back (std::make_pair(row,
+ (*matrix)(row,index)));
+ };
+ };
+};
+
+
+
+template <>
+void
+FilteredMatrix<BlockSparseMatrix<double>,BlockVector<double> >::
+get_column_entries (const unsigned int /*index*/,
+ std::vector<IndexValuePair> &/*column_entries*/,
+ const bool /*matrix_is_symmetric*/) const
+{
+ // presently not implemented, but
+ // should be fairly simple to do
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+template <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::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 <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::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 <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::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 <class Matrix, class Vector>
+typename FilteredMatrix<Matrix,Vector>::value_type
+FilteredMatrix<Matrix,Vector>::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 <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::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 <class Matrix, class Vector>
+typename FilteredMatrix<Matrix,Vector>::value_type
+FilteredMatrix<Matrix,Vector>::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 <class Matrix, class Vector>
+void
+FilteredMatrix<Matrix,Vector>::
+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 <class Matrix, class Vector>
+unsigned int
+FilteredMatrix<Matrix,Vector>::memory_consumption () const
+{
+ return (MemoryConsumption::memory_consumption (matrix) +
+ MemoryConsumption::memory_consumption (constraints) +
+ MemoryConsumption::memory_consumption (tmp_vector));
+};
+
+
+
+template <>
+void
+FilteredMatrix<SparseMatrix<double>,Vector<double> >::
+allocate_tmp_vector ()
+{
+ tmp_mutex.acquire ();
+ tmp_vector.reinit (matrix->n());
+ tmp_mutex.release ();
+};
+
+
+
+template <>
+void
+FilteredMatrix<SparseMatrix<float>,Vector<float> >::
+allocate_tmp_vector ()
+{
+ tmp_mutex.acquire ();
+ tmp_vector.reinit (matrix->n());
+ tmp_mutex.release ();
+};
+
+
+
+template <>
+void
+FilteredMatrix<BlockSparseMatrix<double>,BlockVector<double> >::
+allocate_tmp_vector ()
+{
+ tmp_mutex.acquire ();
+ tmp_vector.reinit (matrix->n());
+ tmp_mutex.release ();
+};
+
+
+
+template <>
+void
+FilteredMatrix<BlockSparseMatrix<float>,BlockVector<float> >::
+allocate_tmp_vector ()
+{
+ tmp_mutex.acquire ();
+ tmp_vector.reinit (matrix->n());
+ tmp_mutex.release ();
+};
+
+
+#endif
--- /dev/null
+//---------------------------- 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 <lac/filtered_matrix.templates.h>
+
+
+template class FilteredMatrix<SparseMatrix<double>,Vector<double> >;
+template class FilteredMatrix<BlockSparseMatrix<double>,BlockVector<double> >;
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
############################################################
--- /dev/null
+//----------------------------------------------------------------------
+// $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 <base/quadrature_lib.h>
+#include <base/function_lib.h>
+#include <lac/vector.h>
+#include <lac/vector_memory.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/filtered_matrix.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <grid/grid_generator.h>
+#include <fe/fe_q.h>
+#include <fe/mapping_q.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <vector>
+#include <fstream>
+#include <string>
+
+
+
+
+void
+solve_filtered (std::map<unsigned int,double> &bv,
+ SparseMatrix<double> &A,
+ Vector<double> &u,
+ Vector<double> &f)
+{
+ FilteredMatrix<SparseMatrix<double> > A1 (A);
+ A1.add_constraints (bv);
+
+ SolverControl control (1000, 1.e-10, false, false);
+ PrimitiveVectorMemory<Vector<double> > mem;
+ SolverCG<Vector<double> > solver (control, mem);
+ PreconditionJacobi<FilteredMatrix<SparseMatrix<double> > > prec;
+ prec.initialize (A1, 1.2);
+
+ Vector<double> f1 (f.size());
+ f1 = f;
+ A1.apply_constraints (f1, true);
+
+ solver.solve (A1, u, f1, prec);
+
+ for (typename std::map<unsigned int,double>::const_iterator i=bv.begin();
+ i!=bv.end(); ++i)
+ Assert (std::fabs(u(i->first) - i->second) < 1e-8,
+ ExcInternalError());
+};
+
+
+
+template <int dim>
+void
+solve_eliminated (std::map<unsigned int,double> &bv,
+ SparseMatrix<double> &A,
+ Vector<double> &u,
+ Vector<double> &f)
+{
+ MatrixTools<dim>::apply_boundary_values (bv, A, u, f);
+
+ SolverControl control (1000, 1.e-10, false, false);
+ PrimitiveVectorMemory<Vector<double> > mem;
+ SolverCG<Vector<double> > solver (control, mem);
+ PreconditionJacobi<> prec;
+ prec.initialize (A, 1.2);
+
+ solver.solve (A, u, f, prec);
+};
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+
+ CosineFunction<dim> cosine;
+
+ if (dim==2)
+ GridGenerator::hyper_ball(tr);
+ else
+ GridGenerator::hyper_cube(tr, -1,1);
+
+ tr.refine_global (5-dim);
+
+ MappingQ<dim> mapping(2);
+ FE_Q<dim> element(1);
+ QGauss4<dim> quadrature;
+
+ DoFHandler<dim> dof(tr);
+ dof.distribute_dofs(element);
+
+ FEValues<dim> fe (mapping, element, quadrature,
+ update_values | update_gradients
+ | update_q_points | update_JxW_values);
+
+ vector <unsigned int> global_dofs (element.dofs_per_cell);
+ vector <double> function (quadrature.n_quadrature_points);
+
+ Vector<double> 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<double> A(A_pattern);
+
+ DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+ const DoFHandler<dim>::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<quadrature.n_quadrature_points;++k)
+ {
+ double dx = fe.JxW (k);
+
+ for (unsigned int i=0;i<element.dofs_per_cell;++i)
+ {
+ const double v = fe.shape_value (i,k);
+ const Tensor<1,dim> grad_v = fe.shape_grad(i,k);
+
+ double rhs = dx * v * (function[k]);
+
+ f(global_dofs[i]) += rhs;
+ for (unsigned int j=0;j<element.dofs_per_cell;++j)
+ {
+ const Tensor<1,dim> 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<unsigned int,double> 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<unsigned int,double>::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<double> 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<double> u_eliminated (dof.n_dofs ());
+ solve_eliminated<dim> (bv, A, u_eliminated, f);
+
+ // output and check
+ for (unsigned int i=0; i<dof.n_dofs(); ++i)
+ {
+ deallog << u_filtered(i) << std::endl;
+ Assert (std::fabs(u_filtered(i) - u_eliminated(i)) < 1e-8,
+ ExcInternalError());
+ };
+}
+
+
+int main ()
+{
+ ofstream logfile ("filtered_matrix.output");
+ logfile.precision (2);
+ logfile.setf(ios::fixed);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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