]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
implement ADI, awaits testing
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 17 Jun 2011 22:26:30 +0000 (22:26 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 17 Jun 2011 22:26:30 +0000 (22:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@23840 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-20/step-20.cc
deal.II/include/deal.II/lac/relaxation_block.h
deal.II/include/deal.II/lac/relaxation_block.templates.h

index 73646d2e33323badf260c61e5c6ba0454de2a92b..c24995dbfcc3edfa686fd1784406f2dcda26fc2c 100644 (file)
@@ -34,6 +34,7 @@
                                 // inverse of a matrix by calling an
                                 // iterative solver.
 #include <deal.II/lac/iterative_inverse.h>
+#include <lac/schur_matrix.h>
 
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
@@ -726,86 +727,6 @@ void MixedLaplaceProblem<dim>::assemble_system ()
                                  // rather only comment on
                                  // implementational aspects.
 
-
-                                 // @sect4{The <code>SchurComplement</code> 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 <code>Subscriptor</code> class and
-                                 // that as mentioned above it stores
-                                 // pointers to the entire block
-                                 // matrix and the inverse of the mass
-                                 // matrix block using
-                                 // <code>SmartPointer</code> objects.
-                                 //
-                                 // The <code>vmult</code> 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 <code>vmult</code>
-                                 // function required these vectors,
-                                 // not the <code>vmult</code> function
-                                 // itself), we allocate them
-                                 // directly, rather than going
-                                 // through the <code>VectorMemory</code>
-                                 // 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 <code>mutable</code>.
-                                 //
-                                 // The rest of the (short)
-                                 // implementation of this class is
-                                 // straightforward if you know the
-                                 // order of matrix-vector
-                                 // multiplications performed by the
-                                 // <code>vmult</code> function:
-class SchurComplement : public Subscriptor
-{
-  public:
-    SchurComplement (const BlockSparseMatrix<double> &A,
-                     const IterativeInverse<Vector<double> > &Minv);
-
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-    const SmartPointer<const IterativeInverse<Vector<double> > > m_inverse;
-    
-    mutable Vector<double> tmp1, tmp2;
-};
-
-
-SchurComplement::SchurComplement (const BlockSparseMatrix<double> &A,
-                                  const IterativeInverse<Vector<double> > &Minv)
-                :
-                system_matrix (&A),
-                m_inverse (&Minv),
-                tmp1 (A.block(0,0).m()),
-                tmp2 (A.block(0,0).m())
-{}
-
-
-void SchurComplement::vmult (Vector<double>       &dst,
-                             const Vector<double> &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 <code>ApproximateSchurComplement</code> class template}
 
                                  // The third component of our solver
@@ -901,7 +822,14 @@ void MixedLaplaceProblem<dim>::solve ()
   m_inverse.solver.select("cg");
   static ReductionControl inner_control(1000, 0., 1.e-13);
   m_inverse.solver.set_control(inner_control);
-  
+
+  SchurComplement<IterativeInverse<Vector<double> >,
+    SparseMatrix<double>,
+    SparseMatrix<double>,
+    SparseMatrix<double> >
+  schur_complement (m_inverse, system_matrix.block(1,0),
+                   system_matrix.block(1,0));
+
   Vector<double> tmp (solution.block(0).size());
 
                                    // Now on to the first
@@ -925,8 +853,6 @@ void MixedLaplaceProblem<dim>::solve ()
     schur_rhs -= system_rhs.block(1);
 
     
-    SchurComplement
-      schur_complement (system_matrix, m_inverse);
     
     ApproximateSchurComplement
       approximate_schur_complement (system_matrix);
index 15726a3cce92f1555b554e9b599583a4d7154ed1..a379da45f3ebfd39ee11f3f5980514ec06979581 100644 (file)
@@ -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:
+                                         *
+                                         * <ol>
+                                         * <li>If the length of the
+                                         * vector is zero, then the
+                                         * relaxation method will be
+                                         * exectued from first to
+                                         * last block.
+                                         * <li> 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.
+                                         * <li> 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.
+                                         * </ol>
+                                         */
+       std::vector<std::vector<unsigned int> > order;
     };
 
                                     /**
index 791fbb2da6433ebcd451c3f6bb4a1bacb0a4dc11..3abd229e0f2c1f439552b4a2abbdcc21b3588285 100644 (file)
@@ -177,29 +177,39 @@ RelaxationBlock<MATRIX,inverse_type>::do_step (
   const MATRIX &M=*this->A;
   Vector<number2> 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;bi<n_blocks;++bi)
+  
+  for (unsigned int perm=0; perm<n_permutations;++perm)
     {
-      const unsigned int block = backward ? (n_blocks - bi - 1) : bi;
-      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_cell<bs; ++row_cell, ++row)
+      for (unsigned int bi=0;bi<n_blocks;++bi)
        {
-         b_cell(row_cell) = src(*row);
-         for (typename MATRIX::const_iterator entry = M.begin(*row);
-              entry != M.end(*row); ++entry)
-           b_cell(row_cell) -= entry->value() * 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_cell<bs; ++row_cell, ++row)
+           {
+             b_cell(row_cell) = src(*row);
+             for (typename MATRIX::const_iterator entry = M.begin(*row);
+                  entry != M.end(*row); ++entry)
+               b_cell(row_cell) -= entry->value() * 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_cell<bs; ++row_cell, ++row)
+           dst(*row) = prev(*row) + additional_data->relaxation * 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_cell<bs; ++row_cell, ++row)
-       dst(*row) = prev(*row) + additional_data->relaxation * x_cell(row_cell);
     }
 }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.