From: Patrick Esser Date: Fri, 16 Sep 2016 09:37:06 +0000 (+0200) Subject: parallelize block inversion X-Git-Tag: v8.5.0-rc1~665^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F3121%2Fhead;p=dealii.git parallelize block inversion --- diff --git a/include/deal.II/lac/relaxation_block.h b/include/deal.II/lac/relaxation_block.h index cbe0eb7208..e8f1f330dd 100644 --- a/include/deal.II/lac/relaxation_block.h +++ b/include/deal.II/lac/relaxation_block.h @@ -249,6 +249,12 @@ protected: * Control information. */ SmartPointer > additional_data; + +private: + /** + * Computes (the inverse of) a range of blocks. + */ + void block_kernel(const size_type block_begin, const size_type block_end); }; diff --git a/include/deal.II/lac/relaxation_block.templates.h b/include/deal.II/lac/relaxation_block.templates.h index 038803cde6..cb6e2f271e 100644 --- a/include/deal.II/lac/relaxation_block.templates.h +++ b/include/deal.II/lac/relaxation_block.templates.h @@ -94,66 +94,74 @@ inline void RelaxationBlock::invert_diagblocks () { - const MatrixType &M=*A; - FullMatrix M_cell; - if (this->same_diagonal()) { Assert(false, ExcNotImplemented()); } else { - for (size_type block=0; blockblock_list.n_rows(); ++block) - { - const size_type bs = additional_data->block_list.row_length(block); - M_cell.reinit(bs, bs); + // compute blocks in parallel + parallel::apply_to_subranges(0, this->additional_data->block_list.n_rows(), + std_cxx11::bind(&RelaxationBlock::block_kernel, this, + std_cxx11::_1, std_cxx11::_2), + 16); + } + this->inverses_computed(true); +} - // Copy rows for this block - // into the matrix for the - // diagonal block - SparsityPattern::iterator row - = additional_data->block_list.begin(block); - for (size_type row_cell=0; row_cellcolumn()); - entry != M.end(row->column()); ++entry) - { - const size_type column = entry->column(); - const size_type col_cell = additional_data->block_list.row_position(block, column); - if (col_cell != numbers::invalid_size_type) - M_cell(row_cell, col_cell) = entry->value(); - } - } - // Now M_cell contains the - // diagonal block. Now - // store it and its - // inverse, if so requested. - if (this->store_diagonals()) - { - this->diagonal(block).reinit(bs, bs); - this->diagonal(block) = M_cell; - } - switch (this->inversion) + +template +inline +void +RelaxationBlock::block_kernel (const size_type block_begin, const size_type block_end) +{ + const MatrixType &M=*(this->A); + FullMatrix M_cell; + + for (size_type block = block_begin; block < block_end; ++block) + { + const size_type bs = this->additional_data->block_list.row_length(block); + M_cell.reinit(bs, bs); + + // Copy rows for this block into the matrix for the diagonal block + SparsityPattern::iterator row + = this->additional_data->block_list.begin(block); + for (size_type row_cell=0; row_cellcolumn()); + entry != M.end(row->column()); ++entry) { - case PreconditionBlockBase::gauss_jordan: - this->inverse(block).reinit(bs, bs); - this->inverse(block).invert(M_cell); - break; - case PreconditionBlockBase::householder: - this->inverse_householder(block).initialize(M_cell); - break; - case PreconditionBlockBase::svd: - this->inverse_svd(block).reinit(bs, bs); - this->inverse_svd(block) = M_cell; - this->inverse_svd(block).compute_inverse_svd(additional_data->threshold); - break; - default: - Assert(false, ExcNotImplemented()); + const size_type column = entry->column(); + const size_type col_cell = this->additional_data->block_list.row_position(block, column); + if (col_cell != numbers::invalid_size_type) + M_cell(row_cell, col_cell) = entry->value(); } } + // Now M_cell contains the diagonal block. Now store it and its + // inverse, if so requested. + if (this->store_diagonals()) + { + this->diagonal(block).reinit(bs, bs); + this->diagonal(block) = M_cell; + } + switch (this->inversion) + { + case PreconditionBlockBase::gauss_jordan: + this->inverse(block).reinit(bs, bs); + this->inverse(block).invert(M_cell); + break; + case PreconditionBlockBase::householder: + this->inverse_householder(block).initialize(M_cell); + break; + case PreconditionBlockBase::svd: + this->inverse_svd(block).reinit(bs, bs); + this->inverse_svd(block) = M_cell; + this->inverse_svd(block).compute_inverse_svd(this->additional_data->threshold); + break; + default: + Assert(false, ExcNotImplemented()); + } } - this->inverses_computed(true); } namespace internal