From: kanschat Date: Wed, 17 Dec 2008 21:26:46 +0000 (+0000) Subject: replace a restriction by a warning in the documentation X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7e0afec6aa5d9e0380848437915827c9a0663d74;p=dealii-svn.git replace a restriction by a warning in the documentation git-svn-id: https://svn.dealii.org/trunk@17981 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/block_matrix_array.h b/deal.II/lac/include/lac/block_matrix_array.h index 0131929cf4..b3eedb75e9 100644 --- a/deal.II/lac/include/lac/block_matrix_array.h +++ b/deal.II/lac/include/lac/block_matrix_array.h @@ -399,10 +399,6 @@ class BlockMatrixArray : public Subscriptor /*@}*/ -/*! @addtogroup Preconditioners - *@{ - */ - /** * Inversion of a block-triangular matrix. @@ -416,9 +412,10 @@ class BlockMatrixArray : public Subscriptor * performed by the vmult() member function. * * @note While block indices may be duplicated (see BlockMatrixArray) - * to add blocks, this is not allowed for diagonal blocks, since + * to add blocks, this has to be used with caution, since * summing up the inverse of two blocks does not yield the inverse of - * the sum. + * the sum. While the latter would be desirable, we can only perform + * the first. * * The implementation may be a little clumsy, but it should be * sufficient as long as the block sizes are much larger than the @@ -456,6 +453,8 @@ class BlockMatrixArray : public Subscriptor * method for solving. * @until Error * + * + * @ingroup Preconditioners * @author Guido Kanschat, 2001, 2005 */ template @@ -614,7 +613,6 @@ class BlockTrianglePrecondition bool backward; }; -/*@}*/ #ifndef DOXYGEN //--------------------------------------------------------------------------- diff --git a/deal.II/lac/source/block_matrix_array.cc b/deal.II/lac/source/block_matrix_array.cc index 0ee615bcf9..4b779c5052 100644 --- a/deal.II/lac/source/block_matrix_array.cc +++ b/deal.II/lac/source/block_matrix_array.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006 by the deal.II authors +// Copyright (C) 2005, 2006, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -307,26 +307,32 @@ BlockTrianglePrecondition::do_row ( m = this->entries.begin(); typename std::vector::Entry>::const_iterator end = this->entries.end(); - typename std::vector::Entry>::const_iterator - diagonal = end; + std::vector::Entry>::const_iterator> + diagonals; Vector* p_aux = this->mem->alloc(); Vector& aux = *p_aux; aux.reinit(dst.block(row_num), true); - + + // Loop over all entries, since + // they are not ordered by rows. for (; m != end ; ++m) { const unsigned int i=m->row; + // Ignore everything not in + // this row if (i != row_num) continue; const unsigned int j=m->col; + // Only use the lower (upper) + // triangle for forward + // (backward) substitution if (((j > i) && !backward) || ((j < i) && backward)) continue; if (j == i) { - Assert (diagonal == end, ExcMultipleDiagonal(j)); - diagonal = m; + diagonals.push_back(m); } else { if (m->transpose) m->matrix->Tvmult(aux, dst.block(j)); @@ -335,13 +341,37 @@ BlockTrianglePrecondition::do_row ( dst.block(i).add (-1 * m->prefix, aux); } } - Assert (diagonal != end, ExcNoDiagonal(row_num)); - - if (diagonal->transpose) - diagonal->matrix->Tvmult(aux, dst.block(row_num)); + Assert (diagonals.size() != 0, ExcNoDiagonal(row_num)); + + // Inverting the diagonal block is + // simple, if there is only one + // matrix + if (diagonals.size() == 1) + { + if (diagonals[0]->transpose) + diagonals[0]->matrix->Tvmult(aux, dst.block(row_num)); + else + diagonals[0]->matrix->vmult(aux, dst.block(row_num)); + dst.block(row_num).equ (diagonals[0]->prefix, aux); + } else - diagonal->matrix->vmult(aux, dst.block(row_num)); - dst.block(row_num).equ (diagonal->prefix, aux); + { + aux = 0.; + for (unsigned int i=0;iprefix); + if (m->transpose) + m->matrix->Tvmult_add(aux, dst.block(row_num)); + else + m->matrix->vmult_add(aux, dst.block(row_num)); + aux.scale(m->prefix); + } + dst.block(row_num) = aux; + } this->mem->free(p_aux); }