From f57a81eb3cea9e52f0bf7868cb987cf2abfc2ebf Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 10 Sep 2008 15:56:31 +0000 Subject: [PATCH] Updated the Trilinos classes so that a parallel use is now also possible for rectangular matrices. Added a copy operation to deal.II vectors. git-svn-id: https://svn.dealii.org/trunk@16788 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/block_vector.h | 42 +++++++++ .../lac/include/lac/block_vector.templates.h | 18 ++++ .../lac/include/lac/trilinos_block_vector.h | 1 + deal.II/lac/include/lac/vector.h | 25 +++++ deal.II/lac/include/lac/vector.templates.h | 24 +++++ deal.II/lac/source/trilinos_sparse_matrix.cc | 94 +++++++++++-------- 6 files changed, 164 insertions(+), 40 deletions(-) diff --git a/deal.II/lac/include/lac/block_vector.h b/deal.II/lac/include/lac/block_vector.h index d948a30941..90958f5ea3 100644 --- a/deal.II/lac/include/lac/block_vector.h +++ b/deal.II/lac/include/lac/block_vector.h @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -26,6 +27,15 @@ DEAL_II_NAMESPACE_OPEN +#ifdef DEAL_II_USE_TRILINOS +namespace TrilinosWrappers +{ + class Vector; +} +#endif + + + /*! @addtogroup Vectors *@{ */ @@ -129,6 +139,17 @@ class BlockVector : public BlockVectorBase > BlockVector (const BlockVector &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 > BlockVector & operator= (const Vector &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 num_blocks blocks of @@ -396,6 +426,7 @@ class BlockVector : public BlockVectorBase > /*----------------------- Inline functions ----------------------------------*/ + template template BlockVector::BlockVector (const std::vector &n, @@ -467,6 +498,17 @@ BlockVector::operator = (const BlockVector &v) } +#ifdef DEAL_II_USE_TRILINOS +template +inline +BlockVector & +BlockVector::operator = (const TrilinosWrappers::BlockVector &v) +{ + BaseClass::operator = (v); + return *this; +} +#endif + template void BlockVector::scale (const value_type factor) { diff --git a/deal.II/lac/include/lac/block_vector.templates.h b/deal.II/lac/include/lac/block_vector.templates.h index 89f82a6f93..0a6e1e8b78 100644 --- a/deal.II/lac/include/lac/block_vector.templates.h +++ b/deal.II/lac/include/lac/block_vector.templates.h @@ -62,6 +62,24 @@ BlockVector::BlockVector (const BlockVector& v) #endif + +#ifdef DEAL_II_USE_TRILINOS + +template +BlockVector::BlockVector (const TrilinosWrappers::BlockVector &v) +{ + this->block_indices = v.get_block_indices(); + this->components.resize(this->n_blocks()); + + for (unsigned int i=0; in_blocks(); ++i) + this->components[i] = v.block(i); + + BaseClass::collect_sizes(); +} + +#endif + + template void BlockVector::reinit (const unsigned int n_bl, const unsigned int bl_sz, diff --git a/deal.II/lac/include/lac/trilinos_block_vector.h b/deal.II/lac/include/lac/trilinos_block_vector.h index b97a4230ca..8674754723 100644 --- a/deal.II/lac/include/lac/trilinos_block_vector.h +++ b/deal.II/lac/include/lac/trilinos_block_vector.h @@ -231,6 +231,7 @@ namespace TrilinosWrappers * Exception */ DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize); + /** * Exception */ diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h index b0f670401a..014610859d 100644 --- a/deal.II/lac/include/lac/vector.h +++ b/deal.II/lac/include/lac/vector.h @@ -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 & + operator = (const TrilinosWrappers::Vector &v); +#endif + /** * Test for equality. This function * assumes that the present vector and diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index fd38939bd3..1dcabdf37d 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -660,6 +660,30 @@ Vector::operator = (const PETScWrappers::MPI::Vector &v) #endif +#ifdef DEAL_II_USE_TRILINOS + +template +Vector & +Vector::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 template bool diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 9488901c6b..50b3cd382f 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -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 (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 (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 (new Epetra_FECrsMatrix(Copy, row_map, (int*)const_cast(&(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 (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 (new Epetra_FECrsMatrix(Copy, row_map, col_map, (int*)const_cast(&(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 values; - std::vector row_indices; + std::vector values; + std::vector row_indices; - for (unsigned int row=0; rowInsertGlobalValues (row, row_length, - &values[0], &row_indices[0]); - } - } + for (unsigned int row=0; rowInsertGlobalValues (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 + // 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 (new Epetra_FECrsMatrix(Copy, row_map, &n_entries_per_row[0], false)); - else - matrix = std::auto_ptr - (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 - (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 values; std::vector 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(&value), &trilinos_j); + if (ierr > 0) + std::cout << ierr << " " << m() << " " << n() << std::endl; AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j)); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } -- 2.39.5