From: kronbichler Date: Fri, 7 Nov 2008 12:31:41 +0000 (+0000) Subject: Updated a few functions in Trilinos sparse matrix and vector classes to make more... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4cda0247e22080c138edc7d661d86cebe1a46075;p=dealii-svn.git Updated a few functions in Trilinos sparse matrix and vector classes to make more of the tests work. git-svn-id: https://svn.dealii.org/trunk@17511 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/trilinos_block_vector.h b/deal.II/lac/include/lac/trilinos_block_vector.h index daee99f3a3..be98f073da 100644 --- a/deal.II/lac/include/lac/trilinos_block_vector.h +++ b/deal.II/lac/include/lac/trilinos_block_vector.h @@ -504,7 +504,7 @@ namespace TrilinosWrappers * elements. */ BlockVector (const BlockVector &V); - + /** * Creates a block vector * consisting of @@ -516,7 +516,40 @@ namespace TrilinosWrappers * reinit of the blocks. */ BlockVector (const unsigned int num_blocks); - + + /** + * Constructor. Set the number of + * blocks to n.size() and + * initialize each block with + * n[i] zero elements. + * + * References BlockVector.reinit(). + */ + BlockVector (const std::vector &N); + + /** + * Constructor. Set the number of + * blocks to + * n.size(). Initialize the + * vector with the elements + * pointed to by the range of + * iterators given as second and + * third argument. Apart from the + * first argument, this + * constructor is in complete + * analogy to the respective + * constructor of the + * std::vector class, but the + * first argument is needed in + * order to know how to subdivide + * the block vector into + * different blocks. + */ + template + BlockVector (const std::vector &n, + const InputIterator first, + const InputIterator end); + /** * Destructor. Clears memory */ @@ -732,6 +765,37 @@ namespace TrilinosWrappers + inline + BlockVector::BlockVector (const std::vector &N) + { + reinit (N); + } + + + + template + BlockVector::BlockVector (const std::vector &n, + const InputIterator first, + const InputIterator end) + { + // first set sizes of blocks, but + // don't initialize them as we will + // copy elements soon + reinit (n, true); + InputIterator start = first; + for (unsigned int b=0; b(n[b])); + + for (unsigned int i=0; iblock(b)(i) = *start; + } + Assert (start == end, ExcIteratorRangeDoesNotMatchVectorSize()); + } + + + inline BlockVector::BlockVector (const unsigned int num_blocks) { diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index b557e10e76..98129b4a00 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -416,30 +416,26 @@ namespace TrilinosWrappers const unsigned int n_max_entries_per_row); /** - * This constructor 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 + * This constructor 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 - * second one specifies the - * number of columns in the - * total matrix. It also provides - * information for the internal - * arrangement in matrix vector - * products, but is not used for - * the distribution of the - * columns – rather, all - * column elements of a row are - * stored on the same - * processor. The vector - * n_entries_per_row specifies - * the number of entries in each - * row of the newly generated + * second one specifies the number of + * columns in the total matrix. It + * also provides information for the + * internal arrangement in matrix + * vector products, but is not used + * for the distribution of the + * columns – rather, all column + * elements of a row are stored on + * the same processor. The vector + * n_entries_per_row + * specifies the number of entries in + * each row of the newly generated * matrix. */ SparseMatrix (const Epetra_Map &InputRowMap, @@ -452,15 +448,29 @@ namespace TrilinosWrappers * #n columns. The resulting matrix * will be completely stored locally. * - * The number of columns entries - * per row is specified as the - * maximum number of entries - * argument. + * The number of columns entries per + * row is specified as the maximum + * number of entries argument. */ SparseMatrix (const unsigned int m, const unsigned int n, const unsigned int n_max_entries_per_row); + /** + * Generate a matrix that is completely + * stored locally, having #m rows and + * #n columns. The resulting matrix + * will be completely stored locally. + * + * The vector + * n_entries_per_row + * specifies the number of entries in + * each row. + */ + SparseMatrix (const unsigned int m, + const unsigned int n, + const std::vector &n_entries_per_row); + /** * Copy constructor. Sets the * calling matrix to be the same diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h index 2f03830f7e..a0babddce0 100755 --- a/deal.II/lac/include/lac/trilinos_vector.h +++ b/deal.II/lac/include/lac/trilinos_vector.h @@ -418,21 +418,28 @@ namespace TrilinosWrappers : map (InputMap) { - vector = std::auto_ptr (new Epetra_FEVector(map, false)); + vector = std::auto_ptr (new Epetra_FEVector(map)); const int min_my_id = map.MinMyGID(); - const unsigned int size = map.NumMyElements(); + const int size = map.NumMyElements(); Assert (map.MaxLID() == size-1, ExcDimensionMismatch(map.MaxLID(), size-1)); + // Need to copy out values, since the + // deal.II might not use doubles, so + // that a direct access is not possible. std::vector indices (size); + std::vector values (size); for (unsigned int i=0; iReplaceGlobalValues (size, &indices[0], + &values[0]); - const int ierr = vector->ReplaceGlobalValues(size, v.begin()+min_my_id, - &indices[0]); AssertThrow (ierr == 0, ExcTrilinosError()); } @@ -637,15 +644,20 @@ namespace TrilinosWrappers map (v.size(), 0, Epetra_SerialComm()) #endif { - vector = std::auto_ptr (new Epetra_FEVector(map, false)); + vector = std::auto_ptr (new Epetra_FEVector(map)); std::vector indices (v.size()); + std::vector values (v.size()); for (unsigned int i=0; iReplaceGlobalValues((int)v.size(), + &indices[0], &values[0]); - const int ierr = vector->ReplaceGlobalValues(v.size(), v.begin(), - &indices[0]); - AssertThrow (ierr == 0, ExcTrilinosError()); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index dfc1be2f26..ea603df8b9 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -167,6 +167,25 @@ namespace TrilinosWrappers int(n_max_entries_per_row), false))) {} + SparseMatrix::SparseMatrix (const unsigned int m, + const unsigned int n, + const std::vector &n_entries_per_row) + : +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + row_map (m, 0, Epetra_MpiComm(MPI_COMM_WORLD)), + col_map (n, 0, Epetra_MpiComm(MPI_COMM_WORLD)), +#else + row_map (m, 0, Epetra_SerialComm()), + col_map (n, 0, Epetra_SerialComm()), +#endif + last_action (Zero), + compressed (true), + matrix (std::auto_ptr + (new Epetra_FECrsMatrix(Copy, row_map, + (int*)const_cast(&(n_entries_per_row[0])), + false))) + {} + SparseMatrix::SparseMatrix (const SparseMatrix &InputMatrix) : Subscriptor(), @@ -609,8 +628,11 @@ 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 the last thing we did was an + // addition, we need to do all the + // data exchange. if (last_action == Add) { int ierr; @@ -619,6 +641,9 @@ namespace TrilinosWrappers AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } + // If our last action was not an + // addition, we tell that at this + // point here. if (last_action != Insert) last_action = Insert; @@ -638,13 +663,19 @@ namespace TrilinosWrappers // command. Otherwise, we're just // able to replace existing entries. if (matrix->Filled() == false) - ierr = matrix->InsertGlobalValues (trilinos_i, 1, - const_cast(&value), - &trilinos_j); + { + ierr = matrix->InsertGlobalValues (trilinos_i, 1, + const_cast(&value), + &trilinos_j); + if (ierr > 0) + ierr = 0; + } else - ierr = matrix->ReplaceGlobalValues (trilinos_i, 1, - const_cast(&value), - &trilinos_j); + { + ierr = matrix->ReplaceGlobalValues (trilinos_i, 1, + const_cast(&value), + &trilinos_j); + } AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j)); AssertThrow (ierr == 0, ExcTrilinosError(ierr));