From: kanschat Date: Tue, 27 Feb 2007 05:52:31 +0000 (+0000) Subject: new permuted block SOR X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1835d389bc5fbbe23b61ff589880359f07535dbb;p=dealii-svn.git new permuted block SOR git-svn-id: https://svn.dealii.org/trunk@14511 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/precondition_block.h b/deal.II/lac/include/lac/precondition_block.h index 99a777f2e7..14b58b3142 100644 --- a/deal.II/lac/include/lac/precondition_block.h +++ b/deal.II/lac/include/lac/precondition_block.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -153,8 +153,8 @@ class PreconditionBlock : public virtual Subscriptor */ void initialize (const MATRIX& A, const AdditionalData parameters); - - +//TODO:[GK] No idea what these are for. Remove if nobody complains + protected: /** * Initialize matrix and block * size for permuted @@ -190,7 +190,7 @@ class PreconditionBlock : public virtual Subscriptor void invert_permuted_diagblocks( const std::vector& permutation, const std::vector& inverse_permutation); - + public: /** * Deletes the inverse diagonal * block matrices if existent, @@ -696,6 +696,18 @@ class PreconditionBlockJacobi : public virtual Subscriptor, * * See PreconditionBlock for requirements on the matrix. * + * Optionally, the entries of the source vector can be treated in the + * order of the indices in the permutation vector set by + * #set_permutation (or the opposite order for Tvmult()). The inverse + * permutation is used for storing elements back into this + * vector. This functionality is automatically enabled after a call to + * set_permutation() with vectors of nonzero size. + * + * @note The diagonal blocks, like the matrix, are not permuted! + * Therefore, the permutation vector can only swap whole blocks. It + * may not change the order inside blocks or swap single indices + * between blocks. + * * @note Instantiations for this template are provided for @ and * @; others can be generated in application programs (see the * section on @ref Instantiations in the manual). @@ -707,6 +719,15 @@ class PreconditionBlockSOR : public virtual Subscriptor, protected PreconditionBlock { public: + /** + * Set the permutation and its + * inverse. These vectors are + * copied into private data, so + * they can be reused or deleted + * after a call to this function. + */ + void set_permutation(const std::vector& permutation, + const std::vector& inverse_permutation); /** * Define number type of matrix. */ @@ -820,9 +841,20 @@ class PreconditionBlockSOR : public virtual Subscriptor, void Tvmult_add (Vector&, const Vector&) const; private: + /** + * The permutation vector + */ + std::vector permutation; + + /** + * The inverse permutation vector + */ + std::vector inverse_permutation; + /** * Actual implementation of the - * preconditioning algorithm. + * preconditioning algorithm + * called by vmult() and vmult_add(). * * The parameter @p adding does * not have any function, yet. @@ -834,7 +866,9 @@ class PreconditionBlockSOR : public virtual Subscriptor, /** * Actual implementation of the - * preconditioning algorithm. + * preconditioning algorithm + * called by Tvmult() and + * Tvmult_add(). * * The parameter @p adding does * not have any function, yet. @@ -844,6 +878,36 @@ class PreconditionBlockSOR : public virtual Subscriptor, const Vector&, const bool adding) const; + /** + * Actual implementation of the + * preconditioning algorithm + * called by vmult() and + * vmult_add() if #permutation + * vectors are not empty. + * + * The parameter @p adding does + * not have any function, yet. + */ + template + void do_permuted_vmult (Vector&, + const Vector&, + const bool adding) const; + + /** + * Actual implementation of the + * preconditioning algorithm + * called by Tvmult() and + * Tvmult_add() if #permutation + * vectors are not empty. + * + * The parameter @p adding does + * not have any function, yet. + */ + template + void do_permuted_Tvmult (Vector&, + const Vector&, + const bool adding) const; + }; diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index c8808876d6..c015278a6b 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -527,6 +527,23 @@ void PreconditionBlockJacobi /*--------------------- PreconditionBlockSOR -----------------------*/ + +template +void PreconditionBlockSOR::set_permutation ( + const std::vector& p, + const std::vector& i) +{ + Assert (p.size() == i.size(), ExcDimensionMismatch(p.size(), i.size())); + permutation.resize(p.size()); + inverse_permutation.resize(p.size()); + for (unsigned int k=0;k template void PreconditionBlockSOR::do_vmult ( @@ -683,7 +700,7 @@ void PreconditionBlockSOR::do_Tvmult ( FullMatrix M_cell(this->blocksize); for (int icell=this->nblocks-1; icell>=0 ; --icell) { - unsigned int cell = (unsigned int) icell; +// unsigned int cell = (unsigned int) icell; // Collect upper triangle for (row = block_start, row_cell = 0; row_cellblocksize; @@ -708,14 +725,16 @@ void PreconditionBlockSOR::do_Tvmult ( } M_cell.householder(b_cell); M_cell.backward(x_cell,b_cell); + + block_end -= this->blocksize; + block_start += this->blocksize; + // distribute x_cell to dst - for (row=cell*this->blocksize, row_cell=0; + for (row=block_end, row_cell=0; row_cellblocksize; ++row_cell, ++row) dst(row)=this->relaxation*x_cell(row_cell); - block_end -= this->blocksize; - block_start += this->blocksize; } } else @@ -739,24 +758,255 @@ void PreconditionBlockSOR::do_Tvmult ( b_cell(row_cell)=b_cell_row; } this->inverse(cell).vmult(x_cell, b_cell); + + block_end -= this->blocksize; + block_start += this->blocksize; + // distribute x_cell to dst - for (row=cell*this->blocksize, row_cell=0; + for (row=block_end, row_cell=0; + row_cellblocksize; + ++row_cell, ++row) + dst(row)=this->relaxation*x_cell(row_cell); + } +} + + +template +template +void PreconditionBlockSOR::do_permuted_vmult ( + Vector &dst, + const Vector &src, + const bool) const +{ + // introduce the following typedef + // since in the use of exceptions, + // strict C++ requires us to + // specify them fully as they are + // from a template dependent base + // class. thus, we'd have to write + // PreconditionBlock::ExcNoMatrixGivenToUse, + // which is lengthy, but also poses + // some problems to the + // preprocessor due to the comma in + // the template arg list. we could + // then wrap the whole thing into + // parentheses, but that creates a + // parse error for gcc for the + // exceptions that do not take + // args... + typedef PreconditionBlock BaseClass; + + Assert (this->A!=0, ExcNotInitialized()); + + const MATRIX &M=*this->A; + Assert (permutation.size() == M.m(), ExcDimensionMismatch(permutation.size(), M.m())); + + Vector b_cell(this->blocksize), x_cell(this->blocksize); + + // cell_row, cell_column are the + // numbering of the blocks (cells). + // row_cell, column_cell are the local + // numbering of the unknowns in the + // blocks. + // row, column are the global numbering + // of the unkowns. + unsigned int row, row_cell, block_start=0; + number2 b_cell_row; + + if (!this->inverses_ready()) + { + FullMatrix M_cell(this->blocksize); + for (unsigned int cell=0; cell < this->nblocks; ++cell) + { + for (row = permutation[block_start], row_cell = 0; + row_cell < this->blocksize; + ++row_cell, ++row) + { + const typename MATRIX::const_iterator row_end = M.end(row); + typename MATRIX::const_iterator entry = M.begin(row); + + b_cell_row=src(row); + for (; entry != row_end; ++entry) + { + const unsigned int column = entry->column(); + if (inverse_permutation[column] < block_start) + b_cell_row -= entry->value() * dst(column); + else if (column < block_start + this->blocksize) + { + const unsigned int column_cell = column - block_start; + M_cell(row_cell, column_cell) = entry->value(); + } + } + b_cell(row_cell)=b_cell_row; + } + M_cell.householder(b_cell); + M_cell.backward(x_cell,b_cell); + // distribute x_cell to dst + for (row=permutation[block_start], row_cell=0; + row_cellblocksize; + ++row_cell, ++row) + dst(row)=this->relaxation*x_cell(row_cell); + + block_start+=this->blocksize; + } + } + else + for (unsigned int cell=0; cell < this->nblocks; ++cell) + { + for (row = permutation[block_start], row_cell = 0; + row_cellblocksize; + ++row_cell, ++row) + { + const typename MATRIX::const_iterator row_end = M.end(row); + typename MATRIX::const_iterator entry = M.begin(row); + + b_cell_row=src(row); + for (; entry != row_end; ++entry) + { + const unsigned int column = entry->column(); + if (inverse_permutation[column] < block_start) + b_cell_row -= entry->value() * dst(column); + } + b_cell(row_cell)=b_cell_row; + } + this->inverse(permutation[block_start]/this->blocksize).vmult(x_cell, b_cell); + // distribute x_cell to dst + for (row=permutation[block_start], row_cell=0; row_cellblocksize; ++row_cell, ++row) dst(row)=this->relaxation*x_cell(row_cell); - block_end -= this->blocksize; - block_start += this->blocksize; + block_start+=this->blocksize; } } + +template +template +void PreconditionBlockSOR::do_permuted_Tvmult ( + Vector &dst, + const Vector &src, + const bool) const +{ + // introduce the following typedef + // since in the use of exceptions, + // strict C++ requires us to + // specify them fully as they are + // from a template dependent base + // class. thus, we'd have to write + // PreconditionBlock::ExcNoMatrixGivenToUse, + // which is lengthy, but also poses + // some problems to the + // preprocessor due to the comma in + // the template arg list. we could + // then wrap the whole thing into + // parentheses, but that creates a + // parse error for gcc for the + // exceptions that do not take + // args... + typedef PreconditionBlock BaseClass; + + Assert (this->A!=0, ExcNotInitialized()); + + const MATRIX &M=*this->A; + Assert (permutation.size() == M.m(), ExcDimensionMismatch(permutation.size(), M.m())); + + Vector b_cell(this->blocksize), x_cell(this->blocksize); + + // cell_row, cell_column are the + // numbering of the blocks (cells). + // row_cell, column_cell are the local + // numbering of the unknowns in the + // blocks. + // row, column are the global numbering + // of the unkowns. + unsigned int row, row_cell; + unsigned int block_end=this->blocksize * this->nblocks; + number2 b_cell_row; + + if (!this->inverses_ready()) + { + FullMatrix M_cell(this->blocksize); + for (int icell=this->nblocks-1; icell>=0 ; --icell) + { + const unsigned int block_start = block_end-this->blocksize; + // Collect upper triangle + for (row = permutation[block_start], row_cell = 0; + row_cellblocksize; + ++row_cell, ++row) + { + const typename MATRIX::const_iterator row_end = M.end(row); + typename MATRIX::const_iterator entry = M.begin(row); + + b_cell_row=src(row); + for (; entry != row_end; ++entry) + { + const unsigned int column = entry->column(); + if (inverse_permutation[column] >= block_end) + b_cell_row -= entry->value() * dst(column); + else if (column >= block_start) + { + const unsigned int column_cell = column - block_start; + M_cell(row_cell, column_cell) = entry->value(); + } + } + b_cell(row_cell)=b_cell_row; + } + M_cell.householder(b_cell); + M_cell.backward(x_cell,b_cell); + + // distribute x_cell to dst + for (row=permutation[block_start], row_cell=0; + row_cellblocksize; + ++row_cell, ++row) + dst(row)=this->relaxation*x_cell(row_cell); + + block_end = block_start; + } + } + else + for (int icell = this->nblocks-1; icell >= 0 ; --icell) + { + const unsigned int block_start = block_end-this->blocksize; + for (row=permutation[block_start], row_cell=0; + row_cellblocksize; + ++row_cell, ++row) + { + const typename MATRIX::const_iterator row_end = M.end(row); + typename MATRIX::const_iterator entry = M.begin(row); + + b_cell_row=src(row); + for (; entry != row_end; ++entry) + { + const unsigned int column = entry->column(); + if (inverse_permutation[column] >= block_end) + b_cell_row -= entry->value() * dst(column); + } + b_cell(row_cell)=b_cell_row; + } + this->inverse(block_start/this->blocksize).vmult(x_cell, b_cell); + + // distribute x_cell to dst + for (row=permutation[block_start], row_cell=0; + row_cellblocksize; + ++row_cell, ++row) + dst(row)=this->relaxation*x_cell(row_cell); + + block_end = block_start; + } +} + + template template void PreconditionBlockSOR ::vmult (Vector &dst, const Vector &src) const { - do_vmult(dst, src, false); + if (permutation.size() == 0) + do_vmult(dst, src, false); + else + do_permuted_vmult(dst, src, false); } @@ -766,7 +1016,10 @@ void PreconditionBlockSOR ::vmult_add (Vector &dst, const Vector &src) const { - do_vmult(dst, src, true); + if (permutation.size() == 0) + do_vmult(dst, src, true); + else + do_permuted_vmult(dst, src, true); } @@ -776,7 +1029,10 @@ void PreconditionBlockSOR ::Tvmult (Vector &dst, const Vector &src) const { - do_Tvmult(dst, src, false); + if (permutation.size() == 0) + do_Tvmult(dst, src, false); + else + do_permuted_Tvmult(dst, src, false); } @@ -786,7 +1042,10 @@ void PreconditionBlockSOR ::Tvmult_add (Vector &dst, const Vector &src) const { - do_Tvmult(dst, src, true); + if (permutation.size() == 0) + do_Tvmult(dst, src, true); + else + do_permuted_Tvmult(dst, src, true); }