#define __deal2__mg_matrix_h
#include <lac/vector.h>
+#include <lac/pointer_matrix.h>
#include <lac/sparse_matrix.h>
#include <multigrid/mg_base.h>
#include <base/mg_level_object.h>
+#include <base/std_cxx1x/shared_ptr.h>
DEAL_II_NAMESPACE_OPEN
/*!@addtogroup mg */
/*@{*/
+namespace MG
+{
+
+/**
+ * Multilevel matrix. This matrix stores an MGLevelObject of
+ * PointerMatrixBase objects. It implementes the interface defined in
+ * MGMatrixBase, so that it can be used as a matrix in Multigrid.
+ *
+ * @author Guido Kanschat
+ * @date 2002, 2010
+ */
+ template <class VECTOR = Vector<double> >
+ class Matrix
+ : public MGMatrixBase<VECTOR>
+ {
+ public:
+ /**
+ * Default constructor for an
+ * empty object.
+ */
+ Matrix();
+
+ /**
+ * Constructor setting up
+ * pointers to the matrices in
+ * <tt>M</tt> by calling initialize().
+ */
+ template <class MATRIX>
+ Matrix(const MGLevelObject<MATRIX>& M);
+
+ /**
+ * Initialize the object such
+ * that the level
+ * multiplication uses the
+ * matrices in <tt>M</tt>
+ */
+ template <class MATRIX>
+ void
+ initialize(const MGLevelObject<MATRIX>& M);
+
+ virtual void vmult (const unsigned int level, VECTOR& dst, const VECTOR& src) const;
+ virtual void vmult_add (const unsigned int level, VECTOR& dst, const VECTOR& src) const;
+ virtual void Tvmult (const unsigned int level, VECTOR& dst, const VECTOR& src) const;
+ virtual void Tvmult_add (const unsigned int level, VECTOR& dst, const VECTOR& src) const;
+
+ /**
+ * Memory used by this object.
+ */
+ unsigned int memory_consumption () const;
+ private:
+ MGLevelObject<std_cxx1x::shared_ptr<PointerMatrixBase<VECTOR> > > matrices;
+ };
+
+}
+
/**
+ * @deprecated Use the much simpler class MG::Matrix instead.
+ *
* Multilevel matrix. This class implements the interface defined by
* MGMatrixBase, using MGLevelObject of an arbitrary
* matrix class.
/*----------------------------------------------------------------------*/
+namespace MG
+{
+ template <class VECTOR>
+ template <class MATRIX>
+ inline void
+ Matrix<VECTOR>::initialize (const MGLevelObject<MATRIX>& p)
+ {
+ matrices.resize(p.min_level(), p.max_level());
+ for (unsigned int level=p.min_level();level <= p.max_level();++level)
+ matrices[level] = std_cxx1x::shared_ptr<PointerMatrixBase<VECTOR> > (new_pointer_matrix_base(p[level], VECTOR()));
+ }
+
+
+ template <class VECTOR>
+ template <class MATRIX>
+ Matrix<VECTOR>::Matrix (const MGLevelObject<MATRIX>& p)
+ {
+ initialize(p);
+ }
+
+
+ template <class VECTOR>
+ Matrix<VECTOR>::Matrix ()
+ {}
+
+
+ template <class VECTOR>
+ void
+ Matrix<VECTOR>::vmult (
+ const unsigned int level,
+ VECTOR& dst,
+ const VECTOR& src) const
+ {
+ matrices[level]->vmult(dst, src);
+ }
+
+
+ template <class VECTOR>
+ void
+ Matrix<VECTOR>::vmult_add (
+ const unsigned int level,
+ VECTOR& dst,
+ const VECTOR& src) const
+ {
+ matrices[level]->vmult_add(dst, src);
+ }
+
+
+ template <class VECTOR>
+ void
+ Matrix<VECTOR>::Tvmult (const unsigned int level,
+ VECTOR& dst,
+ const VECTOR& src) const
+ {
+ matrices[level]->Tvmult(dst, src);
+ }
+
+
+ template <class VECTOR>
+ void
+ Matrix<VECTOR>::Tvmult_add (const unsigned int level,
+ VECTOR& dst,
+ const VECTOR& src) const
+ {
+ matrices[level]->Tvmult_add(dst, src);
+ }
+
+
+ template <class VECTOR>
+ unsigned int
+ Matrix<VECTOR>::memory_consumption () const
+ {
+ return sizeof(*this) + matrices->memory_consumption();
+ }
+}
+
+/*----------------------------------------------------------------------*/
+
template <class MATRIX, class VECTOR>
MGMatrix<MATRIX, VECTOR>::MGMatrix (MGLevelObject<MATRIX>* p)
:
};
+namespace MG
+{
+/**
+ * Smoother using relaxation classes.
+ *
+ * A relaxation class is an object that has two member functions,
+ * @code
+ * void step(VECTOR& x, const VECTOR& d) const;
+ * void Tstep(VECTOR& x, const VECTOR& d) const;
+ * @endcode
+ * performing one step of the smoothing scheme.
+ *
+ * This class performs smoothing on each level. The operation can be
+ * controlled by several parameters. First, the relaxation parameter
+ * @p omega is used in the underlying relaxation method. @p steps is
+ * the number of relaxation steps on the finest level (on all levels
+ * if @p variable is off). If @p variable is @p true, the number of
+ * smoothing steps is doubled on each coarser level. This results in a
+ * method having the complexity of the W-cycle, but saving grid
+ * transfers. This is the method proposed by Bramble at al.
+ *
+ * The option @p symmetric switches on alternating between the
+ * smoother and its transpose in each step as proposed by Bramble.
+ *
+ * @p transpose uses the transposed smoothing operation using
+ * <tt>Tstep</tt> instead of the regular <tt>step</tt> of the
+ * relaxation scheme.
+ *
+ * If you are using block matrices, the second @p initialize function
+ * offers the possibility to extract a single block for smoothing. In
+ * this case, the multigrid method must be used only with the vector
+ * associated to that single block.
+ *
+ * @author Guido Kanschat,
+ * @date 2003, 2009, 2010
+ */
+template<class RELAX, class VECTOR>
+class SmootherRelaxation : public MGSmoother<VECTOR>
+{
+ public:
+ /**
+ * Constructor. Sets memory and
+ * smoothing parameters.
+ */
+ SmootherRelaxation(const unsigned int steps = 1,
+ const bool variable = false,
+ const bool symmetric = false,
+ const bool transpose = false);
+
+ /**
+ * Initialize for matrices. This
+ * function initializes the
+ * smoothing operator with the
+ * same smoother for each level.
+ *
+ * @p additional_data is an
+ * object of type
+ * @p RELAX::AdditionalData and
+ * is handed to the
+ * initialization function of the
+ * relaxation method.
+ */
+ template <class MATRIX2>
+ void initialize (const MGLevelObject<MATRIX2>& matrices,
+ const typename RELAX::AdditionalData & additional_data = typename RELAX::AdditionalData());
+
+ template <class MATRIX2, class DATA>
+ void initialize (const MGLevelObject<MATRIX2>& matrices,
+ const MGLevelObject<DATA>& additional_data);
+ /**
+ * Initialize for matrix
+ * blocks. This function
+ * initializes the smoothing
+ * operator with the same
+ * smoother for each level.
+ *
+ * @p additional_data is an
+ * object of type
+ * @p RELAX::AdditionalData and
+ * is handed to the
+ * initialization function of the
+ * relaxation method.
+ */
+// template <class MATRIX2>
+// void initialize (const MGLevelObject<MatrixBlock<MATRIX2> >& matrices,
+// const typename RELAX::AdditionalData & additional_data = typename RELAX::AdditionalData());
+
+ /**
+ * Empty all vectors.
+ */
+ void clear ();
+
+ /**
+ * The actual smoothing method.
+ */
+ virtual void smooth (const unsigned int level,
+ VECTOR& u,
+ const VECTOR& rhs) const;
+
+ /**
+ * Memory used by this object.
+ */
+ unsigned int memory_consumption () const;
+ private:
+ /**
+ * Object containing relaxation
+ * methods.
+ */
+ MGLevelObject<RELAX> smoothers;
+};
+}
+
/**
* Smoother using relaxation classes.
*
transpose = flag;
}
+//----------------------------------------------------------------------//
+
+namespace MG
+{
+ template <class RELAX, class VECTOR>
+ inline
+ SmootherRelaxation<RELAX, VECTOR>::SmootherRelaxation(
+ const unsigned int steps,
+ const bool variable,
+ const bool symmetric,
+ const bool transpose)
+ : MGSmoother<VECTOR>(steps, variable, symmetric, transpose)
+ {}
+
+
+ template <class RELAX, class VECTOR>
+ inline void
+ SmootherRelaxation<RELAX, VECTOR>::clear ()
+ {
+ smoothers.clear();
+ }
+
+
+ template <class RELAX, class VECTOR>
+ template <class MATRIX2>
+ inline void
+ SmootherRelaxation<RELAX, VECTOR>::initialize (
+ const MGLevelObject<MATRIX2>& m,
+ const typename RELAX::AdditionalData& data)
+ {
+ const unsigned int min = m.get_minlevel();
+ const unsigned int max = m.get_maxlevel();
+
+ smoothers.resize(min, max);
+
+ for (unsigned int i=min;i<=max;++i)
+ smoothers[i].initialize(m[i], data);
+ }
+
+
+ template <class RELAX, class VECTOR>
+ template <class MATRIX2, class DATA>
+ inline void
+ SmootherRelaxation<RELAX, VECTOR>::initialize (
+ const MGLevelObject<MATRIX2>& m,
+ const MGLevelObject<DATA>& data)
+ {
+ const unsigned int min = m.get_minlevel();
+ const unsigned int max = m.get_maxlevel();
+
+ Assert (data.get_minlevel() == min,
+ ExcDimensionMismatch(data.get_minlevel(), min));
+ Assert (data.get_maxlevel() == max,
+ ExcDimensionMismatch(data.get_maxlevel(), max));
+
+ smoothers.resize(min, max);
+
+ for (unsigned int i=min;i<=max;++i)
+ smoothers[i].initialize(m[i], data[i]);
+ }
+
+ template <class RELAX, class VECTOR>
+ inline void
+ SmootherRelaxation<RELAX, VECTOR>::smooth(
+ const unsigned int level,
+ VECTOR& u,
+ const VECTOR& rhs) const
+ {
+ unsigned int maxlevel = smoothers.get_maxlevel();
+ unsigned int steps2 = this->steps;
+
+ if (this->variable)
+ steps2 *= (1<<(maxlevel-level));
+
+ bool T = this->transpose;
+ if (this->symmetric && (steps2 % 2 == 0))
+ T = false;
+ if (this->debug > 0)
+ deallog << 'S' << level << ' ';
+
+ for (unsigned int i=0; i<steps2; ++i)
+ {
+ if (T)
+ smoothers[level].Tstep(u, rhs);
+ else
+ smoothers[level].step(u, rhs);
+ if (this->symmetric)
+ T = !T;
+ }
+ }
+
+
+ template <class RELAX, class VECTOR>
+ inline unsigned int
+ SmootherRelaxation<RELAX, VECTOR>::
+ memory_consumption () const
+ {
+ return sizeof(*this)
+ + smoothers.memory_consumption()
+ + this->mem->memory_consumption();
+ }
+}
+
+
//----------------------------------------------------------------------//
template <class MATRIX, class RELAX, class VECTOR>
VECTOR& u,
const VECTOR& rhs) const
{
- unsigned int maxlevel = matrices.get_maxlevel();
+ unsigned int maxlevel = smoothers.get_maxlevel();
unsigned int steps2 = this->steps;
if (this->variable)