From 337e12ae900d02888112ae9242feb66a2d73b4a9 Mon Sep 17 00:00:00 2001 From: guido Date: Sat, 19 Jul 2003 21:03:00 +0000 Subject: [PATCH] new const_iterator and conjugate_add improved git-svn-id: https://svn.dealii.org/trunk@7884 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/precondition_block.h | 397 +++++++++++++++++- .../lac/precondition_block.templates.h | 59 +-- deal.II/lac/include/lac/sparse_matrix_ez.h | 66 ++- 3 files changed, 471 insertions(+), 51 deletions(-) diff --git a/deal.II/lac/include/lac/precondition_block.h b/deal.II/lac/include/lac/precondition_block.h index 86fd63a540..18eb95603e 100644 --- a/deal.II/lac/include/lac/precondition_block.h +++ b/deal.II/lac/include/lac/precondition_block.h @@ -24,6 +24,9 @@ template class FullMatrix; template class Vector; +template +class PreconditionBlockJacobi; + /*! @addtogroup Preconditioners *@{ */ @@ -197,6 +200,12 @@ class PreconditionBlock : public virtual Subscriptor * yields a preconditioner. */ void set_same_diagonal (); + + /** + * Does the matrix use only one + * diagonal block? + */ + bool same_diagonal () const; /** * Stores the inverse of the @@ -232,6 +241,12 @@ class PreconditionBlock : public virtual Subscriptor */ unsigned int block_size () const; + /** + * The number of blocks of the + * matrix. + */ + unsigned int n_blocks() const; + /** * Determine an estimate for the * memory consumption (in bytes) @@ -324,7 +339,12 @@ class PreconditionBlock : public virtual Subscriptor */ bool store_diagonals; + /** + * Number of blocks. + */ + unsigned int nblocks; private: + /** * Storage of the inverse * matrices of the diagonal @@ -350,7 +370,7 @@ class PreconditionBlock : public virtual Subscriptor * Flag for diagonal compression. * @ref set_same_diagonal() */ - bool same_diagonal; + bool var_same_diagonal; }; @@ -358,7 +378,7 @@ class PreconditionBlock : public virtual Subscriptor /** * Block Jacobi preconditioning. * See @ref{PreconditionBlock} for requirements on the matrix. - * @author Ralf Hartmann, Guido Kanschat, 1999, 2000 + * @author Ralf Hartmann, Guido Kanschat, 1999, 2000, 2003 */ template class PreconditionBlockJacobi : public virtual Subscriptor, @@ -371,6 +391,148 @@ class PreconditionBlockJacobi : public virtual Subscriptor, typedef typename MATRIX::value_type number; public: + /** + * STL conforming iterator. + */ + class const_iterator + { + private: + /** + * Accessor class for iterators + */ + class Accessor + { + public: + /** + * Constructor. Since we use + * accessors only for read + * access, a const matrix + * pointer is sufficient. + */ + Accessor (const PreconditionBlockJacobi *matrix, + const unsigned int row); + + /** + * 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 this matrix entry. + */ + inverse_type value() const; + + protected: + /** + * The matrix accessed. + */ + const PreconditionBlockJacobi* matrix; + + /** + * Save block size here + * for further reference. + */ + unsigned int bs; + + /** + * Current block number. + */ + unsigned int a_block; + + /** + * Current block for extracting data. + * + * This is either + * @p{a_block} or zero, + * depending on + * @p{same_diagonal}. + */ + unsigned int u_block; + + /** + * Iterator inside block. + */ + typename FullMatrix::const_iterator b_iterator; + + /** + * End of current block. + */ + typename FullMatrix::const_iterator b_end; + + /** + * Make enclosing class a + * friend. + */ + friend class const_iterator; + }; + + public: + /** + * Constructor. + */ + const_iterator(const PreconditionBlockJacobi* matrix, + const unsigned int row); + + /** + * 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 @p{==}. + */ + 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; + + private: + /** + * Store an object of the + * accessor class. + */ + Accessor accessor; + }; + /** * Make initialization function * publicly available. @@ -443,6 +605,29 @@ class PreconditionBlockJacobi : public virtual Subscriptor, template void Tvmult_add (Vector&, const Vector&) const; + /** + * STL-like iterator with the + * first entry. + */ + const_iterator begin () const; + + /** + * Final iterator. + */ + const_iterator end () const; + + /** + * STL-like iterator with the + * first entry of row @p{r}. + */ + const_iterator begin (const unsigned int r) const; + + /** + * Final iterator of row @p{r}. + */ + const_iterator end (const unsigned int r) const; + + private: /** * Actual implementation of the @@ -723,6 +908,14 @@ PreconditionBlock::empty () const } +template +inline unsigned int +PreconditionBlock::n_blocks () const +{ + return nblocks; +} + + template inline inverse_type PreconditionBlock::el ( @@ -745,4 +938,204 @@ PreconditionBlock::el ( return B(ib, jb); } +//----------------------------------------------------------------------// + +template +inline +PreconditionBlockJacobi::const_iterator::Accessor:: +Accessor (const PreconditionBlockJacobi* matrix, + const unsigned int row) + : + matrix(matrix), + b_iterator(&matrix->inverse(0), 0, 0), + b_end(&matrix->inverse(0), 0, 0) +{ + bs = matrix->block_size(); + a_block = row / bs; + + // This is the end accessor, which + // does not hava a valid block. + if (a_block == matrix->n_blocks()) + return; + + const unsigned int r = row % bs; + + u_block = (matrix->same_diagonal()) ? 0 : a_block; + + b_iterator = matrix->inverse(u_block).begin(r); + b_end = matrix->inverse(u_block).end(); + + Assert (a_block < matrix->n_blocks(), + ExcIndexRange(a_block, 0, matrix->n_blocks())); +} + + +template +inline +unsigned int +PreconditionBlockJacobi::const_iterator::Accessor::row() const +{ + Assert (a_block < matrix->n_blocks(), + ExcIteratorPastEnd()); + + return bs * a_block + b_iterator->row(); +} + + +template +inline +unsigned int +PreconditionBlockJacobi::const_iterator::Accessor::column() const +{ + Assert (a_block < matrix->n_blocks(), + ExcIteratorPastEnd()); + + return bs * a_block + b_iterator->column(); +} + + +template +inline +inverse_type +PreconditionBlockJacobi::const_iterator::Accessor::value() const +{ + Assert (a_block < matrix->n_blocks(), + ExcIteratorPastEnd()); + + return b_iterator->value(); +} + + +template +inline +PreconditionBlockJacobi::const_iterator:: +const_iterator(const PreconditionBlockJacobi *matrix, + const unsigned int row) + : + accessor(matrix, row) +{} + + +template +inline +typename PreconditionBlockJacobi::const_iterator & +PreconditionBlockJacobi::const_iterator::operator++ () +{ + Assert (*this != accessor.matrix->end(), ExcIteratorPastEnd()); + + ++accessor.b_iterator; + if (accessor.b_iterator == accessor.b_end) + { + ++accessor.a_block; + if (!accessor.matrix->same_diagonal()) + ++accessor.u_block; + + if (accessor.a_block < accessor.matrix->n_blocks()) + { + accessor.b_iterator = accessor.matrix->inverse(accessor.a_block).begin(); + accessor.b_end = accessor.matrix->inverse(accessor.a_block).end(); + } + } + return *this; +} + + +template +inline +const typename PreconditionBlockJacobi::const_iterator::Accessor & +PreconditionBlockJacobi::const_iterator::operator* () const +{ + return accessor; +} + + +template +inline +const typename PreconditionBlockJacobi::const_iterator::Accessor * +PreconditionBlockJacobi::const_iterator::operator-> () const +{ + return &accessor; +} + + +template +inline +bool +PreconditionBlockJacobi::const_iterator:: +operator == (const const_iterator& other) const +{ + if (accessor.a_block == accessor.matrix->n_blocks() && + accessor.a_block == other.accessor.a_block) + return true; + + if (accessor.a_block != other.accessor.a_block) + return false; + + return (accessor.row() == other.accessor.row() && + accessor.column() == other.accessor.column()); +} + + +template +inline +bool +PreconditionBlockJacobi::const_iterator:: +operator != (const const_iterator& other) const +{ + return ! (*this == other); +} + + +template +inline +bool +PreconditionBlockJacobi::const_iterator:: +operator < (const const_iterator& other) const +{ + return (accessor.row() < other.accessor.row() || + (accessor.row() == other.accessor.row() && + accessor.column() < other.accessor.column())); +} + + +template +inline +typename PreconditionBlockJacobi::const_iterator +PreconditionBlockJacobi::begin () const +{ + return const_iterator(this, 0); +} + + +template +inline +typename PreconditionBlockJacobi::const_iterator +PreconditionBlockJacobi::end () const +{ + return const_iterator(this, n_blocks() * block_size()); +} + + +template +inline +typename PreconditionBlockJacobi::const_iterator +PreconditionBlockJacobi::begin ( + const unsigned int r) const +{ + Assert (rm(), ExcIndexRange(r,0,A->m())); + return const_iterator(this, r); +} + + + +template +inline +typename PreconditionBlockJacobi::const_iterator +PreconditionBlockJacobi::end ( + const unsigned int r) const +{ + Assert (rm(), ExcIndexRange(r,0,A->m())); + return const_iterator(this, r+1); +} + #endif diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index f339aca3cb..634bfb86ab 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -28,7 +28,7 @@ PreconditionBlock::PreconditionBlock (): blocksize(0), A(0), store_diagonals(false), - same_diagonal(false) + var_same_diagonal(false) {} @@ -50,7 +50,7 @@ void PreconditionBlock::clear () if (var_diagonal.size()!=0) var_diagonal.erase(var_diagonal.begin(), var_diagonal.end()); blocksize = 0; - same_diagonal = false; + var_same_diagonal = false; A = 0; } @@ -69,7 +69,9 @@ void PreconditionBlock::initialize ( Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m())); blocksize=bsize; relaxation = parameters.relaxation; - same_diagonal = parameters.same_diagonal; + var_same_diagonal = parameters.same_diagonal; + nblocks = A->m()/bsize; + if (parameters.invert_diagonal) invert_diagblocks(); } @@ -79,7 +81,7 @@ template const FullMatrix& PreconditionBlock::inverse(unsigned int i) const { - if (same_diagonal) + if (same_diagonal()) return var_inverse[0]; Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size())); @@ -93,7 +95,7 @@ PreconditionBlock::diagonal(unsigned int i) const { Assert(store_diagonals, ExcDiagonalsNotStored()); - if (same_diagonal) + if (same_diagonal()) return var_diagonal[0]; Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size())); @@ -113,7 +115,15 @@ void PreconditionBlock::set_same_diagonal() { Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); - same_diagonal = true; + var_same_diagonal = true; +} + + +template +bool +PreconditionBlock::same_diagonal() const +{ + return var_same_diagonal; } @@ -134,11 +144,9 @@ void PreconditionBlock::invert_diagblocks() const MATRIX &M=*A; Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); - const unsigned int n_cells = M.m()/blocksize; - FullMatrix M_cell(blocksize); - if (same_diagonal) + if (same_diagonal()) { deallog << "PreconditionBlock uses only one diagonal block" << std::endl; // Invert only the first block @@ -176,7 +184,7 @@ void PreconditionBlock::invert_diagblocks() // set the @p{var_inverse} array to the right // size. we could do it like this: - // var_inverse = vector<>(n_cells,FullMatrix<>()) + // var_inverse = vector<>(nblocks,FullMatrix<>()) // but this would involve copying many // FullMatrix objects. // @@ -185,14 +193,14 @@ void PreconditionBlock::invert_diagblocks() if (store_diagonals) { std::vector > - tmp(n_cells, FullMatrix(blocksize)); + tmp(nblocks, FullMatrix(blocksize)); var_diagonal.swap (tmp); } if (true) { std::vector > - tmp(n_cells, FullMatrix(blocksize)); + tmp(nblocks, FullMatrix(blocksize)); var_inverse.swap (tmp); // make sure the tmp object // goes out of scope as @@ -201,7 +209,7 @@ void PreconditionBlock::invert_diagblocks() M_cell.clear (); - for (unsigned int cell=0; cell Assert(this->A!=0, ExcNotInitialized()); const MATRIX &M=*this->A; - const unsigned int n_cells=M.m()/this->blocksize; Vector b_cell(this->blocksize), x_cell(this->blocksize); @@ -292,7 +299,7 @@ void PreconditionBlockJacobi if (!this->inverses_ready()) { FullMatrix M_cell(this->blocksize); - for (unsigned int cell=0; cellblocksize, row_cell=0; row_cellblocksize; @@ -318,7 +325,7 @@ void PreconditionBlockJacobi } } else - for (unsigned int cell=0; cellblocksize, row_cell=0; row_cellblocksize; @@ -412,7 +419,6 @@ void PreconditionBlockSOR::do_vmult ( Assert (this->A!=0, ExcNotInitialized()); const MATRIX &M=*this->A; - const unsigned int n_cells=M.m()/this->blocksize; Vector b_cell(this->blocksize), x_cell(this->blocksize); @@ -429,7 +435,7 @@ void PreconditionBlockSOR::do_vmult ( if (!this->inverses_ready()) { FullMatrix M_cell(this->blocksize); - for (unsigned int cell=0; cellblocksize; @@ -464,7 +470,7 @@ void PreconditionBlockSOR::do_vmult ( } } else - for (unsigned int cell=0; cellblocksize; @@ -522,7 +528,6 @@ void PreconditionBlockSOR::do_Tvmult ( Assert (this->A!=0, ExcNotInitialized()); const MATRIX &M=*this->A; - const unsigned int n_cells=M.m()/this->blocksize; Vector b_cell(this->blocksize), x_cell(this->blocksize); @@ -534,13 +539,13 @@ void PreconditionBlockSOR::do_Tvmult ( // row, column are the global numbering // of the unkowns. unsigned int row, row_cell, block_start = 0; - unsigned int block_end=this->blocksize *n_cells; + unsigned int block_end=this->blocksize *nblocks; number2 b_cell_row; if (!this->inverses_ready()) { FullMatrix M_cell(this->blocksize); - for (int icell=n_cells-1; icell>=0 ; --icell) + for (int icell=nblocks-1; icell>=0 ; --icell) { unsigned int cell = (unsigned int) icell; // Collect upper triangle @@ -578,7 +583,7 @@ void PreconditionBlockSOR::do_Tvmult ( } } else - for (int icell=n_cells-1; icell >=0 ; --icell) + for (int icell=nblocks-1; icell >=0 ; --icell) { unsigned int cell = (unsigned int) icell; for (row=cell*this->blocksize, row_cell=0; @@ -672,10 +677,8 @@ void PreconditionBlockSSOR::vmult (Vector &d Vector cell_src(this->blocksize); Vector cell_dst(this->blocksize); - const unsigned int n_cells = this->A->m()/this->blocksize; - // Multiply with diagonal blocks - for (unsigned int cell=0; cellblocksize; @@ -704,10 +707,8 @@ void PreconditionBlockSSOR::Tvmult (Vector & Vector cell_src(this->blocksize); Vector cell_dst(this->blocksize); - const unsigned int n_cells = this->A->m()/this->blocksize; - // Multiply with diagonal blocks - for (unsigned int cell=0; cellblocksize; diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.h b/deal.II/lac/include/lac/sparse_matrix_ez.h index a5a4561165..ebd4faa72b 100644 --- a/deal.II/lac/include/lac/sparse_matrix_ez.h +++ b/deal.II/lac/include/lac/sparse_matrix_ez.h @@ -1635,9 +1635,9 @@ SparseMatrixEZ::conjugate_add (const MATRIXA& A, // Determine minimal and // maximal row for a column in // advance. -/* - std::vector::minrow(B.n(), B.m()); - std::vector::maxrow(B.n(), 0); + + std::vector minrow(B.n(), B.m()); + std::vector maxrow(B.n(), 0); while (b1 != b_final) { const unsigned int r = b1->row(); @@ -1645,29 +1645,55 @@ SparseMatrixEZ::conjugate_add (const MATRIXA& A, minrow[b1->column()] = r; if (r > maxrow[b1->column()]) maxrow[b1->column()] = r; + ++b1; } -*/ - - b1 = B.begin(); - while (b1 != b_final) + + typename MATRIXA::const_iterator ai = A.begin(); + const typename MATRIXA::const_iterator ae = A.end(); + + while (ai != ae) { - const unsigned int i = b1->row(); - const unsigned int k = b1->column(); - typename MATRIXB::const_iterator b2 = B.begin(); - while (b2 != b_final) + const typename MATRIXA::value_type a = ai->value(); + // Don't do anything if + // this entry is zero. + if (a == 0.) continue; + + // Now, loop over all rows + // having possibly a + // nonzero entry in column + // ai->row() + b1 = B.begin(minrow[ai->row()]); + const typename MATRIXB::const_iterator + be1 = B.end(maxrow[ai->row()]); + const typename MATRIXB::const_iterator + be2 = B.end(maxrow[ai->column()]); + + while (b1 != be1) { - const unsigned int j = b2->row(); - const unsigned int l = b2->column(); - - const typename MATRIXA::value_type a = A.el(k,l); - - if (a != 0.) + const double b1v = b1->value(); + // We need the product + // of both. If it is + // zero, we can save + // the work + if (b1->column() == ai->row() && (b1v != 0.)) { - add (i, j, a * b1->value() * b2->value()); + const unsigned int i = b1->row(); + + typename MATRIXB::const_iterator + b2 = B.begin(minrow[ai->column()]); + while (b2 != be2) + { + if (b2->column() == ai->column()) + { + const unsigned int j = b2->row(); + add (i, j, a * b1v * b2->value()); + } + ++b2; + } } - ++b2; + ++b1; } - ++b1; + ++ai; } } } -- 2.39.5