From: guido Date: Wed, 25 Aug 2004 14:49:07 +0000 (+0000) Subject: new functions by Erik X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e02aaeb17c47c37334c522a28ac8f256e8ec6465;p=dealii-svn.git new functions by Erik git-svn-id: https://svn.dealii.org/trunk@9576 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 0ba5770d6c..7a5b84bfb4 100644 --- a/deal.II/lac/include/lac/precondition_block.h +++ b/deal.II/lac/include/lac/precondition_block.h @@ -162,6 +162,43 @@ class PreconditionBlock : public virtual Subscriptor void initialize (const MATRIX& A, const AdditionalData parameters); + + /** + * Initialize matrix and block + * size for permuted + * preconditioning. Additionally + * to the parameters of the other + * initalize() function, we hand + * over two index vectors with + * the permutation and its + * inverse. For the meaning of + * these vectors see + * PreconditionBlockSOR. + * + * In a second step, the inverses + * of the diagonal blocks may be + * computed. Make sure you use + * invert_permuted_diagblocks() + * to yield consistent data. + * + * Additionally, a relaxation + * parameter for derived classes + * may be provided. + */ + void initialize (const MATRIX& A, + const std::vector& permutation, + const std::vector& inverse_permutation, + const AdditionalData parameters); + + /** + * Replacement of + * invert_diagblocks() for + * permuted preconditioning. + */ + void invert_permuted_diagblocks( + const std::vector& permutation, + const std::vector& inverse_permutation); + /** * Deletes the inverse diagonal * block matrices if existent, diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index 35eebf04ed..71cb0c0c50 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -76,6 +76,140 @@ void PreconditionBlock::initialize ( invert_diagblocks(); } +template +void PreconditionBlock::initialize ( + const MATRIX &M, + const std::vector& permutation, + const std::vector& inverse_permutation, + const AdditionalData parameters) +{ + + const unsigned int bsize = parameters.block_size; + + clear(); + Assert (M.m() == M.n(), ExcMatrixNotSquare()); + A = &M; + Assert (bsize>0, ExcIndexRange(bsize, 1, M.m())); + Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m())); + blocksize=bsize; + relaxation = parameters.relaxation; + var_same_diagonal = parameters.same_diagonal; + nblocks = A->m()/bsize; + + if (parameters.invert_diagonal) + invert_permuted_diagblocks(permutation, inverse_permutation); + +} + +template +void PreconditionBlock::invert_permuted_diagblocks( + const std::vector& permutation, + const std::vector& inverse_permutation) +{ + Assert (A!=0, ExcNotInitialized()); + Assert (blocksize!=0, ExcNotInitialized()); + + const MATRIX &M=*A; + Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); + + FullMatrix M_cell(blocksize); + + if (same_diagonal()) + { + deallog << "PreconditionBlock uses only one diagonal block" << std::endl; + // Invert only the first block + // This is a copy of the code in the + // 'else' part, stripped of the outer loop + if (store_diagonals) + var_diagonal.resize(1); + var_inverse.resize(1); + var_inverse[0].reinit(blocksize, blocksize); + + for (unsigned int row_cell=0; row_cellcolumn() < blocksize) + M_cell(row_cell, entry->column()) = entry->value(); + ++entry; + } + } + if (store_diagonals) + var_diagonal[0] = M_cell; + var_inverse[0].invert(M_cell); + } + else + { + // 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. + + // set the @p{var_inverse} array to the right + // size. we could do it like this: + // var_inverse = vector<>(nblocks,FullMatrix<>()) + // but this would involve copying many + // FullMatrix objects. + // + // the following is a neat trick which + // avoids copying + if (store_diagonals) + { + std::vector > + tmp(nblocks, FullMatrix(blocksize)); + var_diagonal.swap (tmp); + } + + if (true) + { + std::vector > + tmp(nblocks, FullMatrix(blocksize)); + var_inverse.swap (tmp); + // make sure the tmp object + // goes out of scope as + // soon as possible + }; + + M_cell = 0; + + for (unsigned int cell=0; cellcolumn()column()]column()]-cell_start; + if (column_cell >= blocksize) + continue; + M_cell(row_cell, column_cell) = entry->value(); + } + } + + if (store_diagonals) + var_diagonal[cell] = M_cell; + var_inverse[cell].invert(M_cell); + } + } +} + + template const FullMatrix&