From: kanschat Date: Tue, 27 Feb 2007 20:04:13 +0000 (+0000) Subject: simplify functions and fix bug in PreconditionBlockSSOR X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=22d9f032877697fe76bb85567e423261dacece64;p=dealii-svn.git simplify functions and fix bug in PreconditionBlockSSOR git-svn-id: https://svn.dealii.org/trunk@14513 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 14b58b3142..ad4cffa324 100644 --- a/deal.II/lac/include/lac/precondition_block.h +++ b/deal.II/lac/include/lac/precondition_block.h @@ -840,74 +840,56 @@ class PreconditionBlockSOR : public virtual Subscriptor, template void Tvmult_add (Vector&, const Vector&) const; - private: - /** - * The permutation vector - */ - std::vector permutation; - - /** - * The inverse permutation vector - */ - std::vector inverse_permutation; - + protected: /** - * Actual implementation of the - * preconditioning algorithm - * called by vmult() and vmult_add(). + * Implementation of the forward + * substitution loop called by + * vmult() and vmult_add(). + * + * If a #permutation is set by + * set_permutation(), it will + * automatically be honored by + * this function. * * The parameter @p adding does * not have any function, yet. */ template - void do_vmult (Vector&, - const Vector&, - const bool adding) const; + void forward (Vector&, + const Vector&, + const bool transpose_diagonal, + const bool adding) const; /** - * Actual implementation of the - * preconditioning algorithm - * called by Tvmult() and - * Tvmult_add(). + * Implementation of the backward + * substitution loop called by + * Tvmult() and Tvmult_add(). + * + * If a #permutation is set by + * set_permutation(), it will + * automatically be honored by + * this function. * * The parameter @p adding does * not have any function, yet. */ template - void do_Tvmult (Vector&, - const Vector&, - const bool adding) const; + void backward (Vector&, + const Vector&, + const bool transpose_diagonal, + const bool adding) const; + private: /** - * 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. + * The permutation vector */ - template - void do_permuted_vmult (Vector&, - const Vector&, - const bool adding) const; - + std::vector permutation; + /** - * 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. + * The inverse permutation vector */ - template - void do_permuted_Tvmult (Vector&, - const Vector&, - const bool adding) const; - + std::vector inverse_permutation; + }; diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index c015278a6b..1d901bdcb6 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -546,9 +546,10 @@ void PreconditionBlockSOR::set_permutation ( template template -void PreconditionBlockSOR::do_vmult ( +void PreconditionBlockSOR::forward ( Vector &dst, const Vector &src, + const bool transpose_diagonal, const bool) const { // introduce the following typedef @@ -572,7 +573,12 @@ void PreconditionBlockSOR::do_vmult ( Assert (this->A!=0, ExcNotInitialized()); const MATRIX &M=*this->A; - + const bool permuted = (permutation.size() != 0); + if (permuted) + { + 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 @@ -584,308 +590,74 @@ void PreconditionBlockSOR::do_vmult ( // of the unkowns. unsigned int row, row_cell, block_start=0; number2 b_cell_row; - - if (!this->inverses_ready()) + // The diagonal block if the + // inverses were not precomputed + FullMatrix M_cell(this->blocksize); + + for (unsigned int cell=0; cell < this->nblocks; ++cell) { - FullMatrix M_cell(this->blocksize); - for (unsigned int cell=0; cell < this->nblocks; ++cell) + const unsigned int permuted_block_start = permuted + ? permutation[block_start] + :block_start; + + for (row = permuted_block_start, row_cell = 0; + row_cell < this->blocksize; + ++row_cell, ++row) { - for (row = 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 typename MATRIX::const_iterator row_end = M.end(row); - typename MATRIX::const_iterator entry = M.begin(row); + const unsigned int column = entry->column(); + const unsigned int inverse_permuted_column = permuted + ? inverse_permutation[column] + : column; - b_cell_row=src(row); - for (; entry != row_end; ++entry) + if (inverse_permuted_column < block_start) + b_cell_row -= entry->value() * dst(column); + else if (!this->inverses_ready() && column < block_start + this->blocksize) { - const unsigned int column = entry->column(); - if (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(); - } + const unsigned int column_cell = column - block_start; + if (transpose_diagonal) + M_cell(column_cell, row_cell) = entry->value(); + else + 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=block_start, row_cell=0; - row_cellblocksize; - ++row_cell, ++row) - dst(row)=this->relaxation*x_cell(row_cell); - - block_start+=this->blocksize; + b_cell(row_cell)=b_cell_row; } - } - else - for (unsigned int cell=0; cell < this->nblocks; ++cell) - { - for (row = 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 (column < block_start) - b_cell_row -= entry->value() * dst(column); - } - b_cell(row_cell)=b_cell_row; - } - this->inverse(cell).vmult(x_cell, b_cell); - // distribute x_cell to dst - for (row=cell*this->blocksize, row_cell=0; - row_cellblocksize; - ++row_cell, ++row) - dst(row)=this->relaxation*x_cell(row_cell); - - block_start+=this->blocksize; - } -} - - -template -template -void PreconditionBlockSOR::do_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; - - 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; - 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) + if (this->inverses_ready()) { -// unsigned int cell = (unsigned int) icell; - // Collect upper triangle - for (row = 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 (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); - - block_end -= this->blocksize; - block_start += this->blocksize; - - // distribute x_cell to dst - for (row=block_end, row_cell=0; - row_cellblocksize; - ++row_cell, ++row) - dst(row)=this->relaxation*x_cell(row_cell); - + if (transpose_diagonal) + this->inverse(cell).Tvmult(x_cell, b_cell); + else + this->inverse(cell).vmult(x_cell, b_cell); } - } - else - for (int icell = this->nblocks-1; icell >= 0 ; --icell) - { - unsigned int cell = (unsigned int) icell; - for (row=cell*this->blocksize, 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 (column >= block_end) - b_cell_row -= entry->value() * dst(column); - } - 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=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) + else { - 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; } + + // distribute x_cell to dst + for (row=permuted_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_start+=this->blocksize; - } } template template -void PreconditionBlockSOR::do_permuted_Tvmult ( +void PreconditionBlockSOR::backward ( Vector &dst, const Vector &src, + const bool transpose_diagonal, const bool) const { // introduce the following typedef @@ -909,7 +681,11 @@ void PreconditionBlockSOR::do_permuted_Tvmult ( Assert (this->A!=0, ExcNotInitialized()); const MATRIX &M=*this->A; - Assert (permutation.size() == M.m(), ExcDimensionMismatch(permutation.size(), M.m())); + const bool permuted = (permutation.size() != 0); + if (permuted) + { + Assert (permutation.size() == M.m(), ExcDimensionMismatch(permutation.size(), M.m())); + } Vector b_cell(this->blocksize), x_cell(this->blocksize); @@ -924,76 +700,70 @@ void PreconditionBlockSOR::do_permuted_Tvmult ( unsigned int block_end=this->blocksize * this->nblocks; number2 b_cell_row; - if (!this->inverses_ready()) + FullMatrix M_cell(this->blocksize); + for (unsigned int cell=this->nblocks; cell!=0 ;) { - FullMatrix M_cell(this->blocksize); - for (int icell=this->nblocks-1; icell>=0 ; --icell) + --cell; + const unsigned int block_start = block_end - this->blocksize; + // Collect upper triangle + const unsigned int permuted_block_start = (permutation.size() != 0) + ? permutation[block_start] + :block_start; + for (row = permuted_block_start, row_cell = 0; + row_cellblocksize; + ++row_cell, ++row) { - 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 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(); + const unsigned int inverse_permuted_column = permuted + ? inverse_permutation[column] + : column; + if (inverse_permuted_column >= block_end) + b_cell_row -= entry->value() * dst(column); + else if (!this->inverses_ready() && column >= block_start) { - 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(); - } + const unsigned int column_cell = column - block_start; + // We need the + // transpose of the + // diagonal block, + // so we switch row + // and column + // indices + if (transpose_diagonal) + M_cell(column_cell, row_cell) = entry->value(); + else + M_cell(row_cell, column_cell) = entry->value(); } - b_cell(row_cell)=b_cell_row; } + b_cell(row_cell)=b_cell_row; + } + if (this->inverses_ready()) + { + if (transpose_diagonal) + this->inverse(cell).Tvmult(x_cell, b_cell); + else + this->inverse(cell).vmult(x_cell, b_cell); + } + else + { 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; } + + + // distribute x_cell to dst + for (row=permuted_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; - } } @@ -1003,10 +773,7 @@ void PreconditionBlockSOR ::vmult (Vector &dst, const Vector &src) const { - if (permutation.size() == 0) - do_vmult(dst, src, false); - else - do_permuted_vmult(dst, src, false); + forward(dst, src, false, false); } @@ -1016,10 +783,7 @@ void PreconditionBlockSOR ::vmult_add (Vector &dst, const Vector &src) const { - if (permutation.size() == 0) - do_vmult(dst, src, true); - else - do_permuted_vmult(dst, src, true); + forward(dst, src, false, true); } @@ -1029,10 +793,7 @@ void PreconditionBlockSOR ::Tvmult (Vector &dst, const Vector &src) const { - if (permutation.size() == 0) - do_Tvmult(dst, src, false); - else - do_permuted_Tvmult(dst, src, false); + backward(dst, src, true, false); } @@ -1042,10 +803,7 @@ void PreconditionBlockSOR ::Tvmult_add (Vector &dst, const Vector &src) const { - if (permutation.size() == 0) - do_Tvmult(dst, src, true); - else - do_permuted_Tvmult(dst, src, true); + backward(dst, src, true, true); } @@ -1067,10 +825,11 @@ void PreconditionBlockSSOR::vmult (Vector &d Vector help; help.reinit(dst); - PreconditionBlockSOR::vmult(help, src); + this->forward(help, src, false, false); Vector cell_src(this->blocksize); Vector cell_dst(this->blocksize); + const double scaling = (2.-this->relaxation)/this->relaxation; // Multiply with diagonal blocks for (unsigned int cell=0; cell < this->nblocks; ++cell) @@ -1083,10 +842,10 @@ void PreconditionBlockSSOR::vmult (Vector &d this->diagonal(cell).vmult(cell_dst, cell_src); for (unsigned int row_cell=0; row_cellblocksize; ++row_cell) - help(row+row_cell)=this->relaxation*cell_dst(row_cell); + help(row+row_cell) = scaling * cell_dst(row_cell); } - PreconditionBlockSOR::Tvmult(dst, help); + this->backward(dst, help, false, false); } template @@ -1097,10 +856,11 @@ void PreconditionBlockSSOR::Tvmult (Vector & Vector help; help.reinit(dst); - PreconditionBlockSOR::Tvmult(help, src); + this->backward(help, src, true, false); Vector cell_src(this->blocksize); Vector cell_dst(this->blocksize); + const double scaling = (2.-this->relaxation)/this->relaxation; // Multiply with diagonal blocks for (unsigned int cell=0; cell < this->nblocks; ++cell) @@ -1110,13 +870,13 @@ void PreconditionBlockSSOR::Tvmult (Vector & for (unsigned int row_cell=0; row_cellblocksize; ++row_cell) cell_src(row_cell)=help(row+row_cell); - this->diagonal(cell).vmult(cell_dst, cell_src); + this->diagonal(cell).Tvmult(cell_dst, cell_src); for (unsigned int row_cell=0; row_cellblocksize; ++row_cell) - help(row+row_cell)=this->relaxation*cell_dst(row_cell); + help(row+row_cell) = scaling * cell_dst(row_cell); } - PreconditionBlockSOR::vmult(dst, help); + this->forward(dst, help, true, false); } DEAL_II_NAMESPACE_CLOSE