]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Updated the Trilinos classes so that a parallel use is now also possible for rectangu...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Sep 2008 15:56:31 +0000 (15:56 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Sep 2008 15:56:31 +0000 (15:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@16788 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/block_vector.h
deal.II/lac/include/lac/block_vector.templates.h
deal.II/lac/include/lac/trilinos_block_vector.h
deal.II/lac/include/lac/vector.h
deal.II/lac/include/lac/vector.templates.h
deal.II/lac/source/trilinos_sparse_matrix.cc

index d948a309415be6d2ecf91457f105af533fa643e5..90958f5ea3c711cc5b67679f017160c049dbc632 100644 (file)
@@ -19,6 +19,7 @@
 #include <lac/vector.h>
 #include <lac/block_indices.h>
 #include <lac/block_vector_base.h>
+#include <lac/trilinos_block_vector.h>
 
 #include <cstdio>
 #include <vector>
 DEAL_II_NAMESPACE_OPEN
 
 
+#ifdef DEAL_II_USE_TRILINOS
+namespace TrilinosWrappers
+{
+  class Vector;
+}
+#endif
+
+
+
 /*! @addtogroup Vectors
  *@{
  */
@@ -129,6 +139,17 @@ class BlockVector : public BlockVectorBase<Vector<Number> >
     BlockVector (const BlockVector<OtherNumber> &v);
 #endif    
 
+
+#ifdef DEAL_II_USE_TRILINOS
+                                      /**
+                                       * A copy constructor taking a
+                                       * (parallel) Trilinos block
+                                       * vector and copying it into
+                                       * the deal.II own format.
+                                       */
+    BlockVector (const TrilinosWrappers::BlockVector &v);
+
+#endif
                                      /**
                                       * Constructor. Set the number of
                                       * blocks to
@@ -196,6 +217,15 @@ class BlockVector : public BlockVectorBase<Vector<Number> >
     BlockVector &
     operator= (const Vector<Number> &V);
 
+#ifdef DEAL_II_USE_TRILINOS
+                                     /**
+                                     * A copy constructor from a
+                                     * Trilinos block vector to a
+                                     * deal.II block vector.
+                                     */
+    BlockVector &
+    operator= (const TrilinosWrappers::BlockVector &V);
+#endif
                                      /**
                                      * Reinitialize the BlockVector to
                                      * contain <tt>num_blocks</tt> blocks of
@@ -396,6 +426,7 @@ class BlockVector : public BlockVectorBase<Vector<Number> >
 /*----------------------- Inline functions ----------------------------------*/
 
 
+
 template <typename Number>
 template <typename InputIterator>
 BlockVector<Number>::BlockVector (const std::vector<unsigned int> &n,
@@ -467,6 +498,17 @@ BlockVector<Number>::operator = (const BlockVector<Number2> &v)
 }
 
 
+#ifdef DEAL_II_USE_TRILINOS
+template <typename Number>
+inline
+BlockVector<Number> &
+BlockVector<Number>::operator = (const TrilinosWrappers::BlockVector &v)
+{
+  BaseClass::operator = (v);
+  return *this;
+}
+#endif
+
 template <typename Number>
 void BlockVector<Number>::scale (const value_type factor)
 {
index 89f82a6f93d0374cf64fa32b5120102fa4762941..0a6e1e8b78317d779b78067c0fa312f3044c356f 100644 (file)
@@ -62,6 +62,24 @@ BlockVector<Number>::BlockVector (const BlockVector<OtherNumber>& v)
 
 #endif
 
+
+#ifdef DEAL_II_USE_TRILINOS
+
+template <typename Number>
+BlockVector<Number>::BlockVector (const TrilinosWrappers::BlockVector &v)
+{
+  this->block_indices = v.get_block_indices();
+  this->components.resize(this->n_blocks());
+
+  for (unsigned int i=0; i<this->n_blocks(); ++i)
+    this->components[i] = v.block(i);
+
+  BaseClass::collect_sizes();
+}
+
+#endif
+
+
 template <typename Number>
 void BlockVector<Number>::reinit (const unsigned int n_bl,
                                  const unsigned int bl_sz,
index b97a4230caaf7771edf0a04cdbef82d3eb8458a1..8674754723065d863a5a90b11a80d178dcf277f3 100644 (file)
@@ -231,6 +231,7 @@ namespace TrilinosWrappers
                                           * Exception
                                           */
         DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize);
+
                                          /**
                                           * Exception
                                           */
index b0f670401ad25d19f677f6ad0a5a116c19178a97..014610859d8aa463b390504446784ccb7c76393b 100644 (file)
@@ -400,6 +400,31 @@ class Vector : public Subscriptor
     operator = (const PETScWrappers::MPI::Vector &v);
 #endif
 
+
+#ifdef DEAL_II_USE_TRILINOS
+                                    /**
+                                     * Another copy operator: copy
+                                     * the values from a (sequential
+                                     * or parallel, depending on the
+                                     * underlying compiler) Trilinos
+                                     * wrapper vector class. This
+                                     * operator is only available if
+                                     * Trilinos was detected during
+                                     * configuration time.
+                                      *
+                                      * Note that due to the communication
+                                      * model used in MPI, this operation can
+                                      * only succeed if all processes do it at
+                                      * the same time. I.e., it is not
+                                      * possible for only one process to
+                                      * obtain a copy of a parallel vector
+                                      * while the other jobs do something
+                                      * else.
+                                     */
+    Vector<Number> &
+    operator = (const TrilinosWrappers::Vector &v);
+#endif
+
                                      /**
                                       * Test for equality. This function
                                       * assumes that the present vector and
index fd38939bd3ea11d83873d7ebac6489de3cd5ea09..1dcabdf37db816e5f9c95d012ae281af0bd43cb3 100644 (file)
@@ -660,6 +660,30 @@ Vector<Number>::operator = (const PETScWrappers::MPI::Vector &v)
 #endif
 
 
+#ifdef DEAL_II_USE_TRILINOS
+
+template <typename Number>
+Vector<Number> &
+Vector<Number>::operator = (const TrilinosWrappers::Vector &v)
+{
+  if (v.size() != vec_size)
+    reinit (v.size(), true);
+  if (vec_size != 0)
+    {
+                                       // get a representation of the vector
+                                       // and copy it
+      TrilinosScalar **start_ptr;
+      int ierr = v.vector->ExtractView (&start_ptr);
+      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      
+      std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());
+    }
+
+  return *this;
+}
+
+#endif
+
 template <typename Number>
 template <typename Number2>
 bool
index 9488901c6b1cf88681d6cd708bcbd5afc02f852a..50b3cd382f03e3ac5ab9e5fd595cae14c9823707 100755 (executable)
@@ -85,7 +85,7 @@ namespace TrilinosWrappers
                   row_map (0, 0, Epetra_SerialComm()),
 #endif
                  col_map (row_map),
-                 last_action (Insert),
+                 last_action (Zero),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
                                (new Epetra_FECrsMatrix(Copy, row_map, 0)))
   {}
@@ -95,7 +95,7 @@ namespace TrilinosWrappers
                  :
                   row_map (InputMap),
                  col_map (row_map),
-                 last_action (Insert),
+                 last_action (Zero),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
                                (new Epetra_FECrsMatrix(Copy, row_map, 
                                        int(n_max_entries_per_row), false)))
@@ -106,7 +106,7 @@ namespace TrilinosWrappers
                  :
                   row_map (InputMap),
                  col_map (row_map),
-                 last_action (Insert),
+                 last_action (Zero),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
                    (new Epetra_FECrsMatrix(Copy, row_map, 
                      (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
@@ -119,7 +119,7 @@ namespace TrilinosWrappers
                  :
                   row_map (InputRowMap),
                   col_map (InputColMap),
-                 last_action (Insert),
+                 last_action (Zero),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
                                (new Epetra_FECrsMatrix(Copy, row_map, col_map, 
                                        int(n_max_entries_per_row), false)))
@@ -131,7 +131,7 @@ namespace TrilinosWrappers
                  :
                   row_map (InputRowMap),
                   col_map (InputColMap),
-                 last_action (Insert),
+                 last_action (Zero),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
                    (new Epetra_FECrsMatrix(Copy, row_map, col_map, 
                      (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
@@ -167,11 +167,9 @@ namespace TrilinosWrappers
                                  // sparsity pattern on that
                                  // processor, and only broadcast the
                                  // pattern afterwards.
-    if (row_map.Comm().MyPID() == 0)
-      {
-       Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
-               ExcDimensionMismatch (matrix->NumGlobalRows(),
-                                     sparsity_pattern.n_rows()));
+    Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+           ExcDimensionMismatch (matrix->NumGlobalRows(),
+                                 sparsity_pattern.n_rows()));
     
                                  // Trilinos seems to have a bug for
                                  // rectangular matrices at this point,
@@ -186,24 +184,24 @@ namespace TrilinosWrappers
 //             ExcDimensionMismatch (matrix->NumGlobalCols(),
 //                                   sparsity_pattern.n_cols()));
 
-       std::vector<double> values;
-       std::vector<int>    row_indices;
+    std::vector<double> values;
+    std::vector<int>    row_indices;
     
-       for (unsigned int row=0; row<n_rows; ++row)
-         {
-           const int row_length = sparsity_pattern.row_length(row);
-           row_indices.resize (row_length, -1);
-           values.resize (row_length, 0.);
-           
-           for (int col=0; col < row_length; ++col)
-             row_indices[col] = sparsity_pattern.column_number (row, col);
-           
-           matrix->InsertGlobalValues (row, row_length,
-                                       &values[0], &row_indices[0]);
-         }
-      }
+    for (unsigned int row=0; row<n_rows; ++row)
+      if (row_map.MyGID(row))
+       {
+         const int row_length = sparsity_pattern.row_length(row);
+         row_indices.resize (row_length, -1);
+         values.resize (row_length, 0.);
+
+         for (int col=0; col < row_length; ++col)
+           row_indices[col] = sparsity_pattern.column_number (row, col);
+
+         matrix->InsertGlobalValues (row, row_length,
+                                     &values[0], &row_indices[0]);
+       }
     
-    last_action = Insert;
+    last_action = Zero;
 
                                  // In the end, the matrix needs to
                                  // be compressed in order to be
@@ -250,18 +248,20 @@ namespace TrilinosWrappers
       n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
 
                                  // TODO: There seems to be problem
-                                 // in Epetra when a quadratic matrix
-                                 // is initialized with both row and
+                                 // in Epetra when a matrix is
+                                 // initialized with both row and
                                  // column map. Maybe find something
-                                 // more out about this...
-    if (row_map.SameAs(col_map) == true)
-      matrix = std::auto_ptr<Epetra_FECrsMatrix>
+                                 // more out about this... It could
+                                 // be related to the bug #4123. For
+                                 // the moment, just ignore the
+                                 // column information and generate
+                                 // the matrix as if it were
+                                 // square. The call to
+                                 // GlobalAssemble later will set the
+                                 // correct values.
+    matrix = std::auto_ptr<Epetra_FECrsMatrix>
                (new Epetra_FECrsMatrix(Copy, row_map,
                                        &n_entries_per_row[0], false));
-    else
-      matrix = std::auto_ptr<Epetra_FECrsMatrix>
-               (new Epetra_FECrsMatrix(Copy, row_map, col_map,
-                                       &n_entries_per_row[0], false));
 
     reinit (sparsity_pattern);
   }
@@ -316,9 +316,21 @@ namespace TrilinosWrappers
       n_entries_per_row[(int)row] = 
          dealii_sparse_matrix.get_sparsity_pattern().row_length(row);
 
+                                 // TODO: There seems to be problem
+                                 // in Epetra when a matrix is
+                                 // initialized with both row and
+                                 // column map. Maybe find something
+                                 // more out about this... It could
+                                 // be related to the bug #4123. For
+                                 // the moment, just ignore the
+                                 // column information and generate
+                                 // the matrix as if it were
+                                 // square. The call to
+                                 // GlobalAssemble later will set the
+                                 // correct values.
     matrix = std::auto_ptr<Epetra_FECrsMatrix>
-             (new Epetra_FECrsMatrix(Copy, row_map, col_map,
-                                     &n_entries_per_row[0], true));
+               (new Epetra_FECrsMatrix(Copy, row_map,
+                                       &n_entries_per_row[0], false));
 
     std::vector<double> values;
     std::vector<int> row_indices;
@@ -446,7 +458,7 @@ namespace TrilinosWrappers
 
     Assert (numbers::is_finite(value), 
            ExcMessage("The given value is not finite but either "
-                     "infinite or Not A Number (NaN)"));
+                      "infinite or Not A Number (NaN)"));
 
     if (last_action == Insert)
       {
@@ -455,11 +467,11 @@ namespace TrilinosWrappers
          ierr = matrix->GlobalAssemble (false);
        else
          ierr = matrix->GlobalAssemble(col_map, row_map, false);
-       
+
        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
        last_action = Add;
-      }
+    }
 
                                  // we have to do above actions in any
                                  // case to be consistent with the MPI
@@ -478,6 +490,8 @@ namespace TrilinosWrappers
                                                  const_cast<double*>(&value), 
                                                  &trilinos_j);
 
+    if (ierr > 0)
+      std::cout << ierr << " " << m() << " " << n() << std::endl; 
     AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }

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.