]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added a new function for data exchange to Trilinos vector. Added a block direct solve...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 6 Oct 2008 17:02:37 +0000 (17:02 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 6 Oct 2008 17:02:37 +0000 (17:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@17113 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_block_sparse_matrix.h
deal.II/lac/include/lac/trilinos_block_vector.h
deal.II/lac/include/lac/trilinos_precondition.h
deal.II/lac/source/trilinos_block_vector.cc
deal.II/lac/source/trilinos_precondition.cc
deal.II/lac/source/trilinos_vector.cc

index 4b9d3a761a62a52e15ea4120cd03602cbca006e5..afb5e31cc33b73b37474f7c202103bfe296cbedd 100644 (file)
@@ -116,7 +116,7 @@ namespace TrilinosWrappers
                                         */
       ~BlockSparseMatrix ();
 
-                                       /** 
+                                       /**
                                         * Pseudo copy operator only copying
                                         * empty objects. The sizes of the block
                                         * matrices need to be the same.
index 30ca2aec75ec79dabe8c5b98b1d415650ac971df..daee99f3a3c4c5a2bd3356e68e21f5e139708221 100644 (file)
@@ -39,6 +39,7 @@ namespace TrilinosWrappers
     class BlockVector;
   }
   class BlockVector;
+  class BlockSparseMatrix;
 
 
   namespace MPI
@@ -236,7 +237,50 @@ namespace TrilinosWrappers
                                           */
         void reinit (const unsigned int num_blocks);
 
-                                        /** 
+                                      /**
+                                       * This reinit function is
+                                       * meant to be used for
+                                       * parallel calculations where
+                                       * some non-local data has to
+                                       * be used. The typical
+                                       * situation where one needs
+                                       * this function is the call of
+                                       * the
+                                       * FEValues<dim>::get_function_values
+                                       * function (or of some
+                                       * derivatives) in
+                                       * parallel. Since it is
+                                       * usually faster to retrieve
+                                       * the data in advance, this
+                                       * function can be called
+                                       * before the assembly forks
+                                       * out to the different
+                                       * processors. What this
+                                       * function does is the
+                                       * following: It takes the
+                                       * information in the columns
+                                       * of the given matrix and
+                                       * looks which data couples
+                                       * between the different
+                                       * processors. That data is
+                                       * then queried from the input
+                                       * vector. Note that you should
+                                       * not write to the resulting
+                                       * vector any more, since the
+                                       * some data can be stored
+                                       * several times on different
+                                       * processors, leading to
+                                       * unpredictable results. In
+                                       * particular, such a vector
+                                       * cannot be used for
+                                       * matrix-vector products as
+                                       * for example done during the
+                                       * solution of linear systems.
+                                       */
+       void do_data_exchange (const TrilinosWrappers::BlockSparseMatrix &m,
+                              const BlockVector                         &v);
+
+                                        /**
                                          * Compresses all the components
                                          * after assembling together all
                                          * elements.
@@ -245,7 +289,7 @@ namespace TrilinosWrappers
 
                                         /**
                                          * Returns the state of the
-                                         * matrix, i.e., whether
+                                         * vector, i.e., whether
                                          * compress() needs to be
                                          * called after an operation
                                          * requiring data
index ac8def9bcafeef49e93883427c0f24e4be247506..20932338642f47b2b06b60d14780cfd5ecfba59b 100755 (executable)
@@ -88,7 +88,7 @@ namespace TrilinosWrappers
       friend class SolverBase;
       friend class SolverSaddlePoint;
 
-    protected:
+      //protected:
                                       /**
                                        * This is a pointer to the
                                        * preconditioner object.
@@ -407,8 +407,9 @@ namespace TrilinosWrappers
  * AdditionalData structure for details), and a parameter
  * <tt>overlap</tt> that determines if and how much overlap there
  * should be between the matrix partitions on the various MPI
- * processes. The default settings are 1 for the relaxation parameter,
- * 0 for the diagonal augmentation and 0 for the overlap.
+ * processes.  The default settings are 0 for the additional fill-in, 0
+ * for the absolute augmentation tolerance, 1 for the relative
+ * augmentation tolerance, 0 for the overlap.
  *
  * Note that a parallel applicatoin of the IC preconditioner is
  * actually a block-Jacobi preconditioner with block size equal to the
@@ -533,8 +534,9 @@ namespace TrilinosWrappers
  * AdditionalData structure for details), and a parameter
  * <tt>overlap</tt> that determines if and how much overlap there
  * should be between the matrix partitions on the various MPI
- * processes. The default settings are 1 for the relaxation parameter,
- * 0 for the diagonal augmentation and 0 for the overlap.
+ * processes. The default settings are 0 for the additional fill-in, 0
+ * for the absolute augmentation tolerance, 1 for the relative
+ * augmentation tolerance, 0 for the overlap.
  *
  * Note that a parallel applicatoin of the ILU preconditioner is
  * actually a block-Jacobi preconditioner with block size equal to the
@@ -633,6 +635,60 @@ namespace TrilinosWrappers
       void initialize (const SparseMatrix   &matrix,
                       const AdditionalData &additional_data = AdditionalData());
   };
+
+
+
+/**
+ * A wrapper class for a sparse direct LU decomposition on parallel
+ * blocks for Trilinos matrices. When run in serial, this corresponds
+ * to a direct solve on the matrix.
+ *
+ * The AdditionalData data structure allows to set preconditioner
+ * options.
+ *
+ * Note that a parallel applicatoin of the block direct solve
+ * preconditioner is actually a block-Jacobi preconditioner with block
+ * size equal to the local matrix size. Spoken more technically, this
+ * parallel operation is an <a
+ * href="http://en.wikipedia.org/wiki/Additive_Schwarz_method">additive
+ * Schwarz method</a> with an <em>exact solve</em> as inner solver,
+ * based on the (outer) parallel partitioning.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+  class PreconditionBlockDirect : public PreconditionBase
+  {
+    public:
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional parameters
+                                        * to the preconditioner.
+                                        */      
+      struct AdditionalData
+      {
+                                       /**
+                                       * Constructor.
+                                       */
+       AdditionalData (const unsigned int overlap  = 0);
+       
+                                       /**
+                                       * Block direct parameters and overlap.
+                                       */
+       unsigned int overlap;
+      };
+
+                                       /**
+                                        * Initialize function. Takes
+                                        * the matrix which is used to
+                                        * form the preconditioner, and
+                                        * additional flags if there
+                                        * are any.
+                                        */
+      void initialize (const SparseMatrix   &matrix,
+                      const AdditionalData &additional_data = AdditionalData());
+  };
   
 
 }
index 27f52188f431e93c0f66c78b9b088b403a6bad31..acaa66ffbf5ca2798904e369a8a7db00c683ac22 100644 (file)
@@ -12,6 +12,7 @@
 //---------------------------------------------------------------------------
 
 #include <lac/trilinos_block_vector.h>
+#include <lac/trilinos_block_sparse_matrix.h>
 
 
 #ifdef DEAL_II_USE_TRILINOS
@@ -135,6 +136,29 @@ namespace TrilinosWrappers
 
 
 
+    void
+    BlockVector::do_data_exchange (const TrilinosWrappers::BlockSparseMatrix &m,
+                                  const BlockVector                         &v)
+    {
+      Assert (m.n_block_rows() == v.n_blocks(),
+             ExcDimensionMismatch(m.n_block_rows(),v.n_blocks()));
+      Assert (m.n_block_cols() == v.n_blocks(),
+             ExcDimensionMismatch(m.n_block_cols(),v.n_blocks()));
+
+      if (v.n_blocks() != n_blocks())
+       {
+         block_indices = v.get_block_indices();
+         components.resize(v.n_blocks());
+       }
+
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+       components[i].do_data_exchange(m.block(i,i), v.block(i));
+
+      collect_sizes();
+    }
+
+
+
     void
     BlockVector::compress ()
     {
index 55bc5c1bcd3e37736f32c377dfcb22edd4b05c8e..9b4c4ae03c664fb11cd20179c96d545e864a2945 100755 (executable)
@@ -72,6 +72,7 @@ namespace TrilinosWrappers
     int ierr;
 
     Teuchos::ParameterList parameter_list;
+    parameter_list.set ("relaxation: sweeps", 1);
     parameter_list.set ("relaxation: type", "Jacobi");
     parameter_list.set ("relaxation: damping factor", additional_data.omega);
     parameter_list.set ("relaxation: min diagonal value", 
@@ -118,6 +119,7 @@ namespace TrilinosWrappers
     int ierr;
 
     Teuchos::ParameterList parameter_list;
+    parameter_list.set ("relaxation: sweeps", 1);
     parameter_list.set ("relaxation: type", "symmetric Gauss-Seidel");
     parameter_list.set ("relaxation: damping factor", additional_data.omega);
     parameter_list.set ("relaxation: min diagonal value", 
@@ -165,6 +167,7 @@ namespace TrilinosWrappers
     int ierr;
 
     Teuchos::ParameterList parameter_list;
+    parameter_list.set ("relaxation: sweeps", 1);
     parameter_list.set ("relaxation: type", "Gauss-Seidel");
     parameter_list.set ("relaxation: damping factor", additional_data.omega);
     parameter_list.set ("relaxation: min diagonal value", 
@@ -275,6 +278,44 @@ namespace TrilinosWrappers
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }
 
+
+
+/* ---------------------- PreconditionBlockDirect --------------------- */
+
+  PreconditionBlockDirect::AdditionalData::
+  AdditionalData (const unsigned int overlap)
+                  :
+                  overlap  (overlap)
+  {}
+
+
+
+  void
+  PreconditionBlockDirect::initialize (const SparseMatrix   &matrix,
+                                      const AdditionalData &additional_data)
+  {
+    preconditioner.release();
+
+    preconditioner = Teuchos::rcp (Ifpack().Create ("Amesos", &*matrix.matrix, 
+                                                   additional_data.overlap));
+    Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+                                              "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("schwarz: combine mode", "Add");
+    
+    ierr = preconditioner->SetParameters(parameter_list);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = preconditioner->Initialize();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = preconditioner->Compute();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
 }
 
 DEAL_II_NAMESPACE_CLOSE
index 3501c8818154f14166f61f5a727d9510128ac33a..5e4ad408fee7d0d6d5bf380780bd19be872587bb 100755 (executable)
@@ -215,6 +215,33 @@ namespace TrilinosWrappers
       return *this;
     }
 
+
+
+    void
+    Vector::do_data_exchange (const TrilinosWrappers::SparseMatrix &m,
+                             const Vector                         &v)
+    {
+      Assert (m.matrix->Filled() == true,
+             ExcMessage ("Matrix is not compressed. "
+                         "Cannot find exchange information!"));
+      Assert (v.vector->Map().UniqueGIDs() == true,
+             ExcMessage ("The input vector has overlapping data, "
+                         "which is not allowed."));
+
+      if (vector->Map().SameAs(m.matrix->ColMap()) == false)
+       {
+         map = m.matrix->ColMap();
+         vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
+       }
+
+      Epetra_Import data_exchange (vector->Map(), v.vector->Map());
+      const int ierr = vector->Import(*v.vector, data_exchange, Insert);
+
+      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+      last_action = Insert;
+    }
+
   } /* end of namespace MPI */
 
 

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.