From e79d27737e62b273bca7d6c7fad75f82c51efcef Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 17 Jan 2000 20:21:25 +0000 Subject: [PATCH] Compute inverses at beginning. git-svn-id: https://svn.dealii.org/trunk@2241 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_vanka.h | 58 ++++++-- .../lac/include/lac/sparse_vanka.templates.h | 128 ++++++++++++++---- 2 files changed, 145 insertions(+), 41 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_vanka.h b/deal.II/lac/include/lac/sparse_vanka.h index 1191b6ea0c..718f47fa2a 100644 --- a/deal.II/lac/include/lac/sparse_vanka.h +++ b/deal.II/lac/include/lac/sparse_vanka.h @@ -98,7 +98,7 @@ * whether this might pose some problems in the inversion of the local matrices. * Maybe someone would like to check this. * - * @author Guido Kanschat, documentation mostly by Wolfgang Bangerth; 1999 + * @author Guido Kanschat, documentation and extensions by Wolfgang Bangerth; 1999, 2000 */ template class SparseVanka @@ -114,16 +114,50 @@ class SparseVanka * longer than the Vanka * object. The same is true for * the matrix. + * + * The matrix #M# which is passed + * here may or may not be the + * same matrix for which this + * object shall act as + * preconditioner. In particular, + * it is conceivable that the + * preconditioner is build up for + * one matrix once, but is used + * for subsequent steps in a + * nonlinear process as well, + * where the matrix changes in + * each step slightly. + * + * If #conserve_mem# is #false#, + * then the inverses of the local + * systems are computed now, if + * the flag is #true#, then they + * are computed every time the + * preconditioner is + * applied. This saves some + * memory, but makes + * preconditioning very + * slow. Note also, that if the + * flag is #false#, the the + * contents of the matrix #M# at + * the time of calling this + * constructor are used, while if + * the flag is #true#, then the + * values in #M# at the time of + * preconditioning are used. This + * may lead to different results, + * obviously, of #M# changes. */ SparseVanka(const SparseMatrix &M, - const vector &selected); + const vector &selected, + const bool conserve_memory = false); /** * Destructor. * Delete all allocated matrices. */ ~SparseVanka(); - + /** * Do the preconditioning. * This function takes the residual @@ -134,18 +168,6 @@ class SparseVanka void operator() (Vector &dst, const Vector &src) const; - /** - * Minimize memory consumption. - * Activating this option reduces - * memory needs of the Vanka object - * to nearly zero. You pay for this - * by a high increase of computing - * time, since all local matrices - * are built up and inverted every - * time the Vanka operator is applied. - */ - void conserve_memory(); - /** * Exception */ @@ -175,6 +197,12 @@ class SparseVanka * that are tagged in #selected#. */ mutable vector > > inverses; + + /** + * Compute the inverses of all + * selected diagonal elements. + */ + void compute_inverses (); }; diff --git a/deal.II/lac/include/lac/sparse_vanka.templates.h b/deal.II/lac/include/lac/sparse_vanka.templates.h index 3591447074..8b5ad7ca04 100644 --- a/deal.II/lac/include/lac/sparse_vanka.templates.h +++ b/deal.II/lac/include/lac/sparse_vanka.templates.h @@ -12,15 +12,19 @@ template SparseVanka::SparseVanka(const SparseMatrix &M, - const vector &selected) + const vector &selected, + const bool conserve_mem) : matrix(&M), selected(selected), - conserve_mem(false), + conserve_mem(conserve_mem), inverses(M.m(),0) { Assert (M.m() == M.n(), ExcMatrixNotSquare ()); + + if (conserve_mem == false) + compute_inverses (); } @@ -28,9 +32,9 @@ template SparseVanka::~SparseVanka() { vector > >::iterator i; - for(i=inverses.begin();i!=inverses.end();++i) + for(i=inverses.begin(); i!=inverses.end(); ++i) { - FullMatrix* p = *i; + FullMatrix *p = *i; *i = 0; if (p != 0) delete p; } @@ -38,6 +42,95 @@ SparseVanka::~SparseVanka() +template +void +SparseVanka::compute_inverses () +{ + // first define an alias to the sparsity + // pattern of the matrix, since this + // will be used quite often + const SparseMatrixStruct &structure + = matrix->get_sparsity_pattern(); + + // traverse all rows of the matrix + // which are selected + for (unsigned int row=0; row< matrix->m() ; ++row) + if (selected[row] == true) + { + const unsigned int row_length = structure.row_length(row); + inverses[row] = new FullMatrix (row_length, + row_length); + // mapping between: + // 1 column number of all + // entries in this row, and + // 2 the position within this + // row (as stored in the + // sparsematrixstruct object + // + // since we do not explicitely + // consider nonsysmmetric sparsity + // patterns, the first element + // of each entry simply denotes + // all degrees of freedom that + // couple with #row#. + map local_index; + for (unsigned int i=0; i + (structure.column_number(row, i), i)); + + // Build local matrix and rhs + for (map::const_iterator is=local_index.begin(); + is!=local_index.end(); ++is) + { + // irow loops over all DoFs that + // couple with the present DoF + const unsigned int irow = is->first; + // index of DoF irow in the matrix + // row corresponding to DoF #row#. + // runs between 0 and row_length + const unsigned int i = is->second; + // number of DoFs coupling to + // irow (including irow itself) + const unsigned int irow_length = structure.row_length(irow); + + // copy rhs +// b(i) = src(irow); + + // for all the DoFs that irow + // couples with + for (unsigned int j=0; j::const_iterator js + = local_index.find(col); + // if not, then still use + // this dof to modify the rhs + // + // note that if so, we already + // have copied the entry above + if (js == local_index.end()) +;//b(i) -= matrix->raw_entry(irow,j) * dst(col); + else + // if so, then build the + // matrix out of it + (*inverses[row])(i,js->second) = matrix->raw_entry(irow,j); + }; + }; + + // Compute new values + inverses[row]->gauss_jordan(); + }; +}; + + + template template void @@ -68,24 +161,8 @@ SparseVanka::operator ()(Vector &dst, for (unsigned int row=0; row< matrix->m() ; ++row) if (selected[row] == true) { - // shall we reconstruct the system - // of equations for this DoF? - bool build_matrix = true; const unsigned int row_length = structure.row_length(row); - if (!conserve_mem) - { - if (inverses[row] == 0) - // inverse not yet built - { - FullMatrix *p = new FullMatrix; - inverses[row] = p; - } else { - // inverse already built - build_matrix = false; - } - } - // if we don't store the // inverse matrices, then alias // the entry in the global @@ -95,10 +172,7 @@ SparseVanka::operator ()(Vector &dst, { inverses[row] = &local_matrix; inverses[row]->reinit (row_length, row_length); - } - else - if (build_matrix) - inverses[row]->reinit (row_length, row_length); + }; b.reinit (row_length); x.reinit (row_length); @@ -163,14 +237,16 @@ SparseVanka::operator ()(Vector &dst, else // if so, then build the // matrix out of it - if (build_matrix) + if (conserve_mem == true) (*inverses[row])(i,js->second) = matrix->raw_entry(irow,j); }; }; // Compute new values - if (build_matrix) + if (conserve_mem == true) inverses[row]->gauss_jordan(); + + // apply preconditioner inverses[row]->vmult(x,b); // Distribute new values -- 2.39.5