From: kanschat Date: Tue, 30 Oct 2007 20:23:25 +0000 (+0000) Subject: add scaled matrix and impreove FilteredMatrix X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4465de3221ac19e6c77e283baf8687a89f806981;p=dealii-svn.git add scaled matrix and impreove FilteredMatrix git-svn-id: https://svn.dealii.org/trunk@15397 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/filtered_matrix.h b/deal.II/lac/include/lac/filtered_matrix.h index 0f99834092..17a7dae912 100644 --- a/deal.II/lac/include/lac/filtered_matrix.h +++ b/deal.II/lac/include/lac/filtered_matrix.h @@ -209,6 +209,133 @@ template class FilteredMatrix : public Subscriptor { public: + class const_iterator; + + /** + * Accessor class for iterators + */ + class Accessor + { + /** + * Constructor. Since we use + * accessors only for read + * access, a const matrix + * pointer is sufficient. + */ + Accessor (const FilteredMatrix *matrix, + const unsigned int index); + + public: + /** + * Row number of the element + * represented by this + * object. + */ + unsigned int row() const; + + /** + * Column number of the + * element represented by + * this object. + */ + unsigned int column() const; + + /** + * Value of the right hand + * side for this row. + */ + double value() const; + + private: + /** + * Advance to next entry + */ + void advance (); + + /** + * The matrix accessed. + */ + const FilteredMatrix* matrix; + + /** + * Current row number. + */ + unsigned int index; + /* + * Make enclosing class a + * friend. + */ + friend class const_iterator; + }; + + /** + * STL conforming iterator. + */ + class const_iterator + { + public: + /** + * Constructor. + */ + const_iterator(const FilteredMatrix *matrix, + const unsigned int index); + + /** + * Prefix increment. + */ + const_iterator& operator++ (); + + /** + * Postfix increment. + */ + const_iterator& operator++ (int); + + /** + * Dereferencing operator. + */ + const Accessor& operator* () const; + + /** + * Dereferencing operator. + */ + const Accessor* operator-> () const; + + /** + * Comparison. True, if + * both iterators point to + * the same matrix + * position. + */ + bool operator == (const const_iterator&) const; + /** + * Inverse of ==. + */ + bool operator != (const const_iterator&) const; + + /** + * Comparison operator. Result is + * true if either the first row + * number is smaller or if the row + * numbers are equal and the first + * index is smaller. + */ + bool operator < (const const_iterator&) const; + + /** + * Comparison operator. Compares just + * the other way around than the + * operator above. + */ + bool operator > (const const_iterator&) const; + + private: + /** + * Store an object of the + * accessor class. + */ + Accessor accessor; + }; + /** * Typedef defining a type that * represents a pair of degree of @@ -216,7 +343,11 @@ class FilteredMatrix : public Subscriptor * shall have. */ typedef std::pair IndexValuePair; - + +/** + * @name Constructors and initialization + */ +//@{ /** * Default constructor. You will * have to set the matrix to be @@ -273,6 +404,19 @@ class FilteredMatrix : public Subscriptor template void initialize (const MATRIX &m, bool expect_constrained_source = false); + +//@} +/** + * @name Managing constraints + */ +//@{ + /** + * Add the constraint that the + * value with index i + * should have the value + * v. + */ + void add_constraint (const unsigned int i, const double v); /** * Add a list of constraints to @@ -289,7 +433,7 @@ class FilteredMatrix : public Subscriptor * @p std::vector of * IndexValuePair objects, * but also a - * std::map. + * std::map. * * The second component of these * pairs will only be used in @@ -306,13 +450,17 @@ class FilteredMatrix : public Subscriptor */ template void add_constraints (const ConstraintList &new_constraints); - + /** * Delete the list of constraints * presently in use. */ void clear_constraints (); - +//@} +/** + * Vector operations + */ +//@{ /** * Apply the constraints to a * right hand side vector. This @@ -394,7 +542,23 @@ class FilteredMatrix : public Subscriptor */ void Tvmult_add (VECTOR &dst, const VECTOR &src) const; - +//@} + +/** + * @name Iterators + */ +//@{ + /** + * Iterator to the first + * constraint. + */ + const_iterator begin () const; + /** + * Final iterator. + */ + const_iterator end () const; +//@} + /** * Determine an estimate for the * memory consumption (in bytes) @@ -514,12 +678,154 @@ class FilteredMatrix : public Subscriptor */ void post_filter (const VECTOR &in, VECTOR &out) const; + + friend class Accessor; }; /*@}*/ /*---------------------- Inline functions -----------------------------------*/ +//--------------------------------Iterators--------------------------------------// + +template +inline +FilteredMatrix::Accessor::Accessor( + const FilteredMatrix *matrix, + const unsigned int index) + : + matrix(matrix), + index(index) +{ + Assert (index <= matrix->constraints.size(), + ExcIndexRange(index, 0, matrix->constraints.size())); +} + + + +template +inline +unsigned int +FilteredMatrix::Accessor::row() const +{ + return matrix->constraints[index].first; +} + + + +template +inline +unsigned int +FilteredMatrix::Accessor::column() const +{ + return matrix->constraints[index].first; +} + + + +template +inline +double +FilteredMatrix::Accessor::value() const +{ + return matrix->constraints[index].second; +} + + + +template +inline +void +FilteredMatrix::Accessor::advance() +{ + Assert (index < matrix->constraints.size(), ExcIteratorPastEnd()); + ++index; +} + + + + +template +inline +FilteredMatrix::const_iterator::const_iterator( + const FilteredMatrix *matrix, + const unsigned int index) + : + accessor(matrix, index) +{} + + + +template +inline +typename FilteredMatrix::const_iterator& +FilteredMatrix::const_iterator::operator++ () +{ + accessor.advance(); + return *this; +} + + +template +inline +const typename FilteredMatrix::Accessor & +FilteredMatrix::const_iterator::operator* () const +{ + return accessor; +} + + +template +inline +const typename FilteredMatrix::Accessor * +FilteredMatrix::const_iterator::operator-> () const +{ + return &accessor; +} + + +template +inline +bool +FilteredMatrix::const_iterator:: +operator == (const const_iterator& other) const +{ + return (accessor.index == other.accessor.index + && accessor.matrix == other.accessor.matrix); +} + + +template +inline +bool +FilteredMatrix::const_iterator:: +operator != (const const_iterator& other) const +{ + return ! (*this == other); +} + + + +//------------------------------- FilteredMatrix ---------------------------------------// + +template +inline +typename FilteredMatrix::const_iterator +FilteredMatrix::begin () const +{ + return const_iterator(this, 0); +} + + +template +inline +typename FilteredMatrix::const_iterator +FilteredMatrix::end () const +{ + return const_iterator(this, constraints.size()); +} + + template inline bool @@ -589,6 +895,18 @@ FilteredMatrix::operator = (const FilteredMatrix &fm) +template +inline +void +FilteredMatrix:: +add_constraint (const unsigned int index, const double value) +{ + // add new constraint to end + constraints.push_back(IndexValuePair(index, value)); +} + + + template template inline @@ -670,10 +988,10 @@ FilteredMatrix::post_filter (const VECTOR &in, { // iterate over all constraints and // set value correctly - const_index_value_iterator i = constraints.begin(); + 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); + out(i->first) = in(i->first); } @@ -697,9 +1015,12 @@ FilteredMatrix::vmult (VECTOR& dst, const VECTOR& src) const tmp_mutex.release (); } else - matrix->vmult (dst, src); - // finally do post-filtering - post_filter (src, dst); + { + matrix->vmult (dst, src); + } + + // finally do post-filtering + post_filter (src, dst); } @@ -723,7 +1044,10 @@ FilteredMatrix::Tvmult (VECTOR& dst, const VECTOR& src) const tmp_mutex.release (); } else - matrix->Tvmult (dst, src); + { + matrix->Tvmult (dst, src); + } + // finally do post-filtering post_filter (src, dst); } @@ -749,9 +1073,12 @@ FilteredMatrix::vmult_add (VECTOR& dst, const VECTOR& src) const tmp_mutex.release (); } else - matrix->vmult_add (dst, src); - // finally do post-filtering - post_filter (src, dst); + { + matrix->vmult_add (dst, src); + } + + // finally do post-filtering + post_filter (src, dst); } @@ -775,7 +1102,10 @@ FilteredMatrix::Tvmult_add (VECTOR& dst, const VECTOR& src) const tmp_mutex.release (); } else - matrix->Tvmult_add (dst, src); + { + matrix->Tvmult_add (dst, src); + } + // finally do post-filtering post_filter (src, dst); } diff --git a/deal.II/lac/include/lac/matrix_lib.h b/deal.II/lac/include/lac/matrix_lib.h index a06883415c..9fa2e4bac0 100644 --- a/deal.II/lac/include/lac/matrix_lib.h +++ b/deal.II/lac/include/lac/matrix_lib.h @@ -146,6 +146,69 @@ class ProductMatrix : public PointerMatrixBase }; +/** + * A matrix that is the scaled version of another matrix. + * + * Matrix-vector products of this matrix are composed of those of the + * original matrix and scaling by a constant factor. + * + * @author Guido Kanschat, 2007 + */ +template +class ScaledMatrix : public Subscriptor +{ + public: + /** + * Constructor leaving an + * uninitialized object. + */ + ScaledMatrix (); + /** + * Constructor with initialization. + */ + template + ScaledMatrix (const MATRIX& M, const double factor); + + /** + * Destructor + */ + ~ScaledMatrix (); + /** + * Initialize for use with a new + * matrix and factor. + */ + template + void initialize (const MATRIX& M, const double factor); + + /** + * Delete internal matrix pointer. + */ + void clear (); + + /** + * Matrix-vector product. + */ + void vmult (VECTOR& w, const VECTOR& v) const; + + /** + * Tranposed matrix-vector + * product. + */ + void Tvmult (VECTOR& w, const VECTOR& v) const; + + private: + /** + * The matrix. + */ + PointerMatrixBase* m; + /** + * The scaling factor; + */ + double factor; +}; + + + /** * Poor man's matrix product of two sparse matrices. Stores two * matrices #m1 and #m2 of arbitrary type SparseMatrix and implements @@ -438,6 +501,81 @@ class InverseMatrixRichardson : public Subscriptor /*@}*/ //--------------------------------------------------------------------------- + +template +inline +ScaledMatrix::ScaledMatrix() + : + m(0) +{} + + + +template +template +inline +ScaledMatrix::ScaledMatrix(const MATRIX& mat, const double factor) + : + m(new_pointer_matrix_base(mat, VECTOR())), + factor(factor) +{} + + + +template +template +inline +void +ScaledMatrix::initialize(const MATRIX& mat, const double f) +{ + if (m) delete m; + m = new_pointer_matrix_base(mat, VECTOR()); + factor = f; +} + + + +template +inline +void +ScaledMatrix::clear() +{ + if (m) delete m; + m = 0; +} + + + +template +inline +ScaledMatrix::~ScaledMatrix() +{ + clear (); +} + + +template +inline +void +ScaledMatrix::vmult (VECTOR& w, const VECTOR& v) const +{ + m->vmult(w, v); + w *= factor; +} + + +template +inline +void +ScaledMatrix::Tvmult (VECTOR& w, const VECTOR& v) const +{ + m->Tvmult(w, v); + w *= factor; +} + + +//--------------------------------------------------------------------------- + template ProductMatrix::ProductMatrix () : m1(0), m2(0), mem(0, typeid(*this).name())