From: kanschat Date: Fri, 17 Jun 2011 22:26:30 +0000 (+0000) Subject: implement ADI, awaits testing X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=94972a79134971924e9036584bde8fd4263b78d0;p=dealii-svn.git implement ADI, awaits testing git-svn-id: https://svn.dealii.org/trunk@23840 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc index 73646d2e33..c24995dbfc 100644 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@ -34,6 +34,7 @@ // inverse of a matrix by calling an // iterative solver. #include +#include #include #include @@ -726,86 +727,6 @@ void MixedLaplaceProblem::assemble_system () // rather only comment on // implementational aspects. - - // @sect4{The SchurComplement class template} - - // The next class is the Schur - // complement class. Its rationale - // has also been discussed in length - // in the introduction. The only - // things we would like to note is - // that the class, too, is derived - // from the Subscriptor class and - // that as mentioned above it stores - // pointers to the entire block - // matrix and the inverse of the mass - // matrix block using - // SmartPointer objects. - // - // The vmult function requires - // two temporary vectors that we do - // not want to re-allocate and free - // every time we call this - // function. Since here, we have full - // control over the use of these - // vectors (unlike above, where a - // class called by the vmult - // function required these vectors, - // not the vmult function - // itself), we allocate them - // directly, rather than going - // through the VectorMemory - // mechanism. However, again, these - // member variables do not carry any - // state between successive calls to - // the member functions of this class - // (i.e., we never care what values - // they were set to the last time a - // member function was called), we - // mark these vectors as mutable. - // - // The rest of the (short) - // implementation of this class is - // straightforward if you know the - // order of matrix-vector - // multiplications performed by the - // vmult function: -class SchurComplement : public Subscriptor -{ - public: - SchurComplement (const BlockSparseMatrix &A, - const IterativeInverse > &Minv); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer > system_matrix; - const SmartPointer > > m_inverse; - - mutable Vector tmp1, tmp2; -}; - - -SchurComplement::SchurComplement (const BlockSparseMatrix &A, - const IterativeInverse > &Minv) - : - system_matrix (&A), - m_inverse (&Minv), - tmp1 (A.block(0,0).m()), - tmp2 (A.block(0,0).m()) -{} - - -void SchurComplement::vmult (Vector &dst, - const Vector &src) const -{ - system_matrix->block(0,1).vmult (tmp1, src); - m_inverse->vmult (tmp2, tmp1); - system_matrix->block(1,0).vmult (dst, tmp2); -} - - // @sect4{The ApproximateSchurComplement class template} // The third component of our solver @@ -901,7 +822,14 @@ void MixedLaplaceProblem::solve () m_inverse.solver.select("cg"); static ReductionControl inner_control(1000, 0., 1.e-13); m_inverse.solver.set_control(inner_control); - + + SchurComplement >, + SparseMatrix, + SparseMatrix, + SparseMatrix > + schur_complement (m_inverse, system_matrix.block(1,0), + system_matrix.block(1,0)); + Vector tmp (solution.block(0).size()); // Now on to the first @@ -925,8 +853,6 @@ void MixedLaplaceProblem::solve () schur_rhs -= system_rhs.block(1); - SchurComplement - schur_complement (system_matrix, m_inverse); ApproximateSchurComplement approximate_schur_complement (system_matrix); diff --git a/deal.II/include/deal.II/lac/relaxation_block.h b/deal.II/include/deal.II/lac/relaxation_block.h index 15726a3cce..a379da45f3 100644 --- a/deal.II/include/deal.II/lac/relaxation_block.h +++ b/deal.II/include/deal.II/lac/relaxation_block.h @@ -135,6 +135,41 @@ class RelaxationBlock : * call to LAPACKFullMatrix::compute_inverse_svd(). */ double threshold; + + /** + * The order in which blocks + * should be traversed. This + * vector can initiate + * several modes of + * execution: + * + *
    + *
  1. If the length of the + * vector is zero, then the + * relaxation method will be + * exectued from first to + * last block. + *
  2. If the length is one, + * then the inner vector must + * have the same size as + * the number of blocks. The + * relaxation method is + * applied in the order given + * in this vector. + *
  3. If the outer vector + * has length greater one, + * then the relaxation method + * is applied several times, + * each time in the order + * given by the inner vector + * of the corresponding + * index. This mode can for + * instance be used for ADI + * methods and similar + * direction sweeps. + *
+ */ + std::vector > order; }; /** diff --git a/deal.II/include/deal.II/lac/relaxation_block.templates.h b/deal.II/include/deal.II/lac/relaxation_block.templates.h index 791fbb2da6..3abd229e0f 100644 --- a/deal.II/include/deal.II/lac/relaxation_block.templates.h +++ b/deal.II/include/deal.II/lac/relaxation_block.templates.h @@ -177,29 +177,39 @@ RelaxationBlock::do_step ( const MATRIX &M=*this->A; Vector b_cell, x_cell; + const bool permutation_empty = additional_data->order.size() == 0; + const unsigned int n_permutations = (permutation_empty) + ? 1U : additional_data->order.size(); const unsigned int n_blocks = additional_data->block_list.size(); - for (unsigned int bi=0;biblock_list.block_size(block); - - b_cell.reinit(bs); - x_cell.reinit(bs); - // Collect off-diagonal parts - BlockList::const_iterator row = additional_data->block_list.begin(block); - for (unsigned int row_cell=0; row_cellvalue() * prev(entry->column()); + unsigned int block = backward ? (n_blocks - bi - 1) : bi; + if (!permutation_empty) + block = additional_data->order[perm][block]; + + const unsigned int bs = additional_data->block_list.block_size(block); + + b_cell.reinit(bs); + x_cell.reinit(bs); + // Collect off-diagonal parts + BlockList::const_iterator row = additional_data->block_list.begin(block); + for (unsigned int row_cell=0; row_cellvalue() * prev(entry->column()); + } + // Apply inverse diagonal + this->inverse_vmult(block, x_cell, b_cell); + // Store in result vector + row=additional_data->block_list.begin(block); + for (unsigned int row_cell=0; row_cellrelaxation * x_cell(row_cell); } - // Apply inverse diagonal - this->inverse_vmult(block, x_cell, b_cell); - // Store in result vector - row=additional_data->block_list.begin(block); - for (unsigned int row_cell=0; row_cellrelaxation * x_cell(row_cell); } }