From fb9ea57759c5e5371140a1f7ae918d8dcaa06ce2 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Thu, 21 Aug 2008 18:03:55 +0000 Subject: [PATCH] Extended the rudimentary Trilinos wrappers a little bit. Rectangular matrices do not yet seem to work properly. Need to do some more work on that. git-svn-id: https://svn.dealii.org/trunk@16629 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/trilinos_sparse_matrix.h | 120 ++++++++++++--- deal.II/lac/include/lac/trilinos_vector.h | 8 +- deal.II/lac/source/trilinos_sparse_matrix.cc | 137 ++++++++++++++---- deal.II/lac/source/trilinos_vector.cc | 2 +- 4 files changed, 209 insertions(+), 58 deletions(-) diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index dee47baf38..becabaa033 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -337,29 +337,52 @@ namespace TrilinosWrappers const std::vector &n_entries_per_row); /** - * Destructor. Made virtual so that one - * can use pointers to this class. + * This function is similar to the + * one above, but it now takes two + * different Epetra maps for rows + * and columns. This interface is + * meant to be used for generating + * rectangular matrices, where one + * map specifies the parallel + * distribution of rows and the other + * the one of the columns. This is + * in contrast to the first + * constructor, where the same map + * is used for both the number of + * rows and the number of columns. + * The number + * of columns per row is specified + * by the maximum number of entries. */ - virtual ~SparseMatrix (); + SparseMatrix (const Epetra_Map &InputRowMap, + const Epetra_Map &InputColMap, + const unsigned int n_max_entries_per_row); + /** - * This operator assigns a scalar to a - * matrix. Since this does usually not - * make much sense (should we set all - * matrix entries to this value? Only - * the nonzero entries of the sparsity - * pattern?), this operation is only - * allowed if the actual value to be - * assigned is zero. This operator only - * exists to allow for the obvious - * notation matrix=0, which - * sets all elements of the matrix to - * zero, but keeps the sparsity pattern - * previously used. + * This function is similar to the + * one above, but it now takes two + * different Epetra maps for rows + * and columns. This interface is + * meant to be used for generating + * rectangular matrices, where one + * map specifies the parallel + * distribution of rows and the other + * the one of the columns. The + * vector n_entries_per_row specifies + * the number of entries in each + * row of the newly generated matrix. */ - SparseMatrix & - operator = (const double d); + SparseMatrix (const Epetra_Map &InputRowMap, + const Epetra_Map &InputColMap, + const std::vector &n_entries_per_row); /** + * Destructor. Made virtual so that one + * can use pointers to this class. + */ + virtual ~SparseMatrix (); + + /** * This function initializes the * Trilinos matrix by attaching all * the elements to the sparsity @@ -394,6 +417,19 @@ namespace TrilinosWrappers * e.g. when the grid is refined. */ void reinit (const Epetra_Map &input_map, + const CompressedSparsityPattern &sparsity_pattern); + + /** + * This function is similar to the + * other initialization function above, + * but now also reassigns the matrix + * rows and columns according to + * two user-supplied Epetra maps. + * To be used e.g. for rectangular + * matrices after remeshing. + */ + void reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, const CompressedSparsityPattern &sparsity_pattern); /** @@ -404,6 +440,24 @@ namespace TrilinosWrappers */ void clear (); + /** + * This operator assigns a scalar to a + * matrix. Since this does usually not + * make much sense (should we set all + * matrix entries to this value? Only + * the nonzero entries of the sparsity + * pattern?), this operation is only + * allowed if the actual value to be + * assigned is zero. This operator only + * exists to allow for the obvious + * notation matrix=0, which + * sets all elements of the matrix to + * zero, but keeps the sparsity pattern + * previously used. + */ + SparseMatrix & + operator = (const double d); + /** * Set the element (i,j) * to @p value. @@ -929,7 +983,7 @@ namespace TrilinosWrappers /** * Exception */ - DeclException4 (ExcAccessToNonlocalElement, + DeclException4 (ExcAccessToNonLocalElement, int, int, int, int, << "You tried to access element (" << arg1 << "/" << arg2 << ")" @@ -938,11 +992,33 @@ namespace TrilinosWrappers << " are stored locally and can be accessed."); /** - * The Epetra Trilinos mapping that + * Exception + */ + DeclException2 (ExcAccessToNonPresentElement, + int, int, + << "You tried to access element (" << arg1 + << "/" << arg2 << ")" + << " of a sparse matrix, but it appears to not " + << " exist in the Trilinos sparsity pattern."); + + /** + * Epetra Trilinos mapping of the + * matrix rows that + * assigns parts of the matrix to + * the individual processes. + * TODO: is it possible to only use + * a pointer? + */ + Epetra_Map row_map; + + /** + * Pointer to the user-supplied + * Epetra Trilinos mapping of the + * matrix columns that * assigns parts of the matrix to * the individual processes. */ - Epetra_Map map; + Epetra_Map col_map; /** * A sparse matrix object in @@ -1040,7 +1116,7 @@ namespace TrilinosWrappers inline const_iterator:: - const_iterator(const SparseMatrix *matrix, + const_iterator(const SparseMatrix *matrix, const unsigned int row, const unsigned int index) : diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h index 7998b13a13..39a68a3f5a 100755 --- a/deal.II/lac/include/lac/trilinos_vector.h +++ b/deal.II/lac/include/lac/trilinos_vector.h @@ -162,7 +162,7 @@ namespace TrilinosWrappers /** * Exception */ - DeclException3 (ExcAccessToNonlocalElement, + DeclException3 (ExcAccessToNonLocalElement, int, int, int, << "You tried to access element " << arg1 << " of a distributed vector, but only elements " @@ -824,12 +824,6 @@ namespace TrilinosWrappers * the communicator and data distribution * object common to all * Trilinos objects used by deal.II. - * TODO: we probably only need a pointer - * to the map, since the information - * is provided from outside and there - * is no need to copy the map. - * Especially not when we have many - * vectors based on the same map. */ Epetra_Map map; diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 5ce2e21e09..51d0d58d60 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -80,22 +80,23 @@ namespace TrilinosWrappers SparseMatrix::SparseMatrix () : #ifdef DEAL_II_COMPILER_SUPPORTS_MPI - map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)), + row_map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)), #else - map (0,0,Epetra_SerialComm()), + row_map (0, 0, Epetra_SerialComm()), #endif + col_map (row_map), matrix (std::auto_ptr - (new Epetra_FECrsMatrix(Copy, map, 0))), + (new Epetra_FECrsMatrix(Copy, row_map, 0))), last_action (Insert) {} - SparseMatrix::SparseMatrix (const Epetra_Map &InputMap, const unsigned int n_max_entries_per_row) : - map (InputMap), + row_map (InputMap), + col_map (row_map), matrix (std::auto_ptr - (new Epetra_FECrsMatrix(Copy, map, + (new Epetra_FECrsMatrix(Copy, row_map, int(n_max_entries_per_row), false))), last_action (Insert) {} @@ -103,15 +104,42 @@ namespace TrilinosWrappers SparseMatrix::SparseMatrix (const Epetra_Map &InputMap, const std::vector &n_entries_per_row) : - map (InputMap), + row_map (InputMap), + col_map (row_map), + matrix (std::auto_ptr + (new Epetra_FECrsMatrix(Copy, row_map, + (int*)const_cast(&(n_entries_per_row[0])), + false))), + last_action (Insert) + {} + + SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap, + const Epetra_Map &InputColMap, + const unsigned int n_max_entries_per_row) + : + row_map (InputRowMap), + col_map (InputColMap), + matrix (std::auto_ptr + (new Epetra_FECrsMatrix(Copy, row_map, col_map, + int(n_max_entries_per_row), false))), + last_action (Insert) + {} + + SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap, + const Epetra_Map &InputColMap, + const std::vector &n_entries_per_row) + : + row_map (InputRowMap), + col_map (InputColMap), matrix (std::auto_ptr - (new Epetra_FECrsMatrix(Copy, map, + (new Epetra_FECrsMatrix(Copy, row_map, col_map, (int*)const_cast(&(n_entries_per_row[0])), - true))), + false))), last_action (Insert) {} + SparseMatrix::~SparseMatrix () { } @@ -124,14 +152,14 @@ namespace TrilinosWrappers { unsigned int n_rows = sparsity_pattern.n_rows(); - + Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), ExcDimensionMismatch (matrix->NumGlobalRows(), sparsity_pattern.n_rows())); - + std::vector values(n_max_entries_per_row, 0.); std::vector row_indices(n_max_entries_per_row); - + for (unsigned int row=0; row n_entries_per_row(n_rows); + + for (unsigned int row=0; row (new Epetra_FECrsMatrix ( + Copy, row_map, &n_entries_per_row[0], true)); + + const unsigned int n_max_entries_per_row = *std::max_element ( + &n_entries_per_row[0], &n_entries_per_row[n_rows-1]); + + reinit (sparsity_pattern, n_max_entries_per_row); + } + + + + void + SparseMatrix::reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const CompressedSparsityPattern &sparsity_pattern) + { + + unsigned int n_rows = sparsity_pattern.n_rows(); + matrix.reset(); + row_map = input_row_map; + col_map = input_col_map; + + Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(), + ExcDimensionMismatch (input_row_map.NumGlobalElements(), + sparsity_pattern.n_rows())); + Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(), + ExcDimensionMismatch (input_col_map.NumGlobalElements(), + sparsity_pattern.n_cols())); std::vector n_entries_per_row(n_rows); @@ -198,7 +264,8 @@ namespace TrilinosWrappers n_entries_per_row[(int)row] = sparsity_pattern.row_length(row); matrix = std::auto_ptr - (new Epetra_FECrsMatrix(Copy, map, &n_entries_per_row[0], true)); + (new Epetra_FECrsMatrix(Copy, row_map, col_map, + &n_entries_per_row[0], true)); const unsigned int n_max_entries_per_row = *std::max_element ( &n_entries_per_row[0], &n_entries_per_row[n_rows-1]); @@ -216,13 +283,13 @@ namespace TrilinosWrappers // generate an empty matrix. matrix.reset(); #ifdef DEAL_II_COMPILER_SUPPORTS_MPI - map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)), + row_map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)), #else - map = Epetra_Map (0,0,Epetra_SerialComm()), + row_map = Epetra_Map (0,0,Epetra_SerialComm()), #endif matrix = std::auto_ptr - (new Epetra_FECrsMatrix(Copy, map, 0)); + (new Epetra_FECrsMatrix(Copy, row_map, 0)); } @@ -232,7 +299,11 @@ namespace TrilinosWrappers { // flush buffers int ierr; - ierr = matrix->GlobalAssemble (true); + if (row_map.SameAs(col_map)) + ierr = matrix->GlobalAssemble (); + else + ierr = matrix->GlobalAssemble (row_map, col_map); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); ierr = matrix->OptimizeStorage (); @@ -269,7 +340,11 @@ namespace TrilinosWrappers if (last_action == Add) { int ierr; - ierr = matrix->GlobalAssemble(false); + if (row_map.SameAs(col_map)) + ierr = matrix->GlobalAssemble (false); + else + ierr = matrix->GlobalAssemble(row_map, col_map, false); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); last_action = Insert; } @@ -281,6 +356,7 @@ namespace TrilinosWrappers const_cast(&value), &trilinos_j); + AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j)); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -299,8 +375,12 @@ namespace TrilinosWrappers if (last_action == Insert) { int ierr; - ierr = matrix->GlobalAssemble(false); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + if (row_map.SameAs(col_map)) + ierr = matrix->GlobalAssemble (false); + else + ierr = matrix->GlobalAssemble(row_map, col_map, false); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); last_action = Add; } @@ -322,6 +402,7 @@ namespace TrilinosWrappers const_cast(&value), &trilinos_j); + AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j)); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -341,7 +422,7 @@ namespace TrilinosWrappers // continue. if ((trilinos_i == -1 ) || (trilinos_j == -1)) { - Assert (false, ExcAccessToNonlocalElement(i, j, local_range().first, + Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first, local_range().second)); } else diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc index 18f8afc3a8..8634dbb269 100755 --- a/deal.II/lac/source/trilinos_vector.cc +++ b/deal.II/lac/source/trilinos_vector.cc @@ -37,7 +37,7 @@ namespace TrilinosWrappers AssertThrow ((static_cast(index) >= vector.map.MinMyGID()) && (static_cast(index) <= vector.map.MaxMyGID()), - ExcAccessToNonlocalElement (index, vector.map.MinMyGID(), + ExcAccessToNonLocalElement (index, vector.map.MinMyGID(), vector.map.MaxMyGID()-1)); return *(*(vector.vector))[index]; } -- 2.39.5