From cff99deda5d90c4fa903b9b499727b2d976a1970 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 27 Aug 2008 11:58:05 +0000 Subject: [PATCH] Restructured the internal procedures in Trilinos wrapper classes to avoid to have the same code at two places. git-svn-id: https://svn.dealii.org/trunk@16672 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/trilinos_block_sparse_matrix.h | 16 + .../lac/include/lac/trilinos_sparse_matrix.h | 122 ++++--- .../source/trilinos_block_sparse_matrix.cc | 51 ++- deal.II/lac/source/trilinos_sparse_matrix.cc | 335 ++++++++---------- 4 files changed, 281 insertions(+), 243 deletions(-) diff --git a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h index b2ee84f753..cde0a7a122 100644 --- a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -181,6 +182,21 @@ namespace TrilinosWrappers void reinit (const std::vector &input_maps, const BlockSparsityPattern &block_sparsity_pattern); + /** + * This function initializes the + * Trilinos matrix using the deal.II + * sparse matrix and the entries stored + * therein. It uses a threshold + * to copy only elements whose + * modulus is larger than the + * threshold (so zeros in the + * deal.II matrix can be filtered + * away). + */ + void reinit (const std::vector &input_maps, + const ::dealii::BlockSparseMatrix &deal_ii_sparse_matrix, + const double drop_tolerance=1e-13); + /** * This function calls the compress() * command of all matrices after diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 8d3bdfd5d4..bce08ab59c 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -264,22 +264,23 @@ namespace TrilinosWrappers * This class implements a wrapper to use the Trilinos distributed sparse matrix * class Epetra_FECrsMatrix. This is precisely the kind of matrix we deal with * all the time - we most likely get it from some assembly process, where also - * entries not locally owned might need to written and hence need to be - * forwarded to the owner. This class is designed to be used in a distributed + * entries not locally owned might need to be written and hence need to be + * forwarded to the owner process. + * This class is designed to be used in a distributed * memory architecture with an MPI compiler on the bottom, but works equally * well also for serial processes. The only requirement for this class to - * work is that Trilinos is installed with the respective compiler as a - * basis. + * work is that Trilinos has been installed with the same compiler as is used + * for generating deal.II. * * The interface of this class is modeled after the existing * SparseMatrix class in deal.II. It has almost the same member * functions, and is often exchangable. However, since Trilinos only supports a - * single scalar type (double), it is - * not templated, and only works with doubles. + * single scalar type (double), it is not templated, and only works with + * doubles. * * Note that Trilinos only guarantees that operations do what you expect if the - * functions @p GlobalAssemble has been called - * after matrix assembly. Therefore, you need to call + * functions @p GlobalAssemble has been called after matrix assembly. + * Therefore, you need to call * SparseMatrix::compress() before you actually use the matrix. This also * calls @p FillComplete that compresses the storage format for sparse * matrices by discarding unused elements. Trilinos allows to continue with @@ -336,7 +337,16 @@ namespace TrilinosWrappers /** * Constructor using an Epetra_Map * and a maximum number of nonzero - * matrix entries. + * matrix entries. Note that this + * number does not need to be exact, + * and it is even allowed that + * the actual matrix structure + * has more nonzero entries than + * specified in the constructor. + * However it is still advantageous to + * provide good estimates here since + * this will considerably increase + * the performance of the matrix. */ SparseMatrix (const Epetra_Map &InputMap, const unsigned int n_max_entries_per_row); @@ -403,53 +413,33 @@ namespace TrilinosWrappers * can use pointers to this class. */ virtual ~SparseMatrix (); - - /** - * This function initializes the - * Trilinos matrix by attaching all - * the elements to the sparsity - * pattern provided as deal argument. - * This function uses a user- - * provided maximum number of - * elements per row. If that is - * not directly available, use one of the - * other reinit functions. - */ - void reinit (const SparsityPattern &sparsity_pattern, - const unsigned int n_max_entries_per_row); /** * This function initializes the - * Trilinos matrix by attaching all - * the elements to the sparsity - * pattern provided as deal argument, - * now calculating the maximum number - * of nonzeros from the sparsity - * pattern internally. + * Trilinos matrix with a deal.II + * sparsity pattern, i.e. it + * makes the Trilinos Epetra + * matrix know the position of + * nonzero entries according to + * the sparsity pattern. Note that, + * when using this function, the + * matrix must already be initialized + * with a suitable Epetra_Map that + * describes the distribution of + * the matrix among the MPI + * processes. Otherwise, an + * error will be thrown. */ void reinit (const SparsityPattern &sparsity_pattern); - /** - * This function initializes the - * Trilinos matrix using the deal.II - * sparse matrix and the entries stored - * therein. It uses a threshold - * to copy only elements whose - * modulus is larger than the - * threshold (so zeros in the - * deal.II matrix can be filtered - * away). - */ - void reinit (const Epetra_Map &input_map, - const ::dealii::SparseMatrix &deal_ii_sparse_matrix, - const double drop_tolerance=1e-13); - /** - * This function is similar to the - * other initialization function above, - * but now also reassigns the matrix - * rows according to a user-supplied - * Epetra map. This might be used + * This function is initializes the + * Trilinos Epetra matrix according + * to the specified sparsity_pattern, + * and also reassigns the matrix + * rows to different processes + * according to a user-supplied + * Epetra map. This might be useful * when the matrix structure changes, * e.g. when the grid is refined. */ @@ -462,13 +452,43 @@ namespace TrilinosWrappers * 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. + * To be used for rectangular + * matrices. */ void reinit (const Epetra_Map &input_row_map, const Epetra_Map &input_col_map, const SparsityPattern &sparsity_pattern); + /** + * This function initializes the + * Trilinos matrix using the deal.II + * sparse matrix and the entries stored + * therein. It uses a threshold + * to copy only elements with + * modulus larger than the + * threshold (so zeros in the + * deal.II matrix can be filtered + * away). + */ + void reinit (const Epetra_Map &input_map, + const ::dealii::SparseMatrix &dealii_sparse_matrix, + const double drop_tolerance=1e-13); + + /** + * This function is similar to the + * other initialization function with + * deal.II sparse matrix input above, + * but now takes Epetra maps for both + * the rows and the columns of the + * matrix. + * To be used for rectangular + * matrices. + */ + void reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const ::dealii::SparseMatrix &dealii_sparse_matrix, + const double drop_tolerance=1e-13); + /** * Release all memory and return * to a state just like after diff --git a/deal.II/lac/source/trilinos_block_sparse_matrix.cc b/deal.II/lac/source/trilinos_block_sparse_matrix.cc index 08b3009734..04f43fe0f5 100644 --- a/deal.II/lac/source/trilinos_block_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_block_sparse_matrix.cc @@ -22,8 +22,8 @@ namespace TrilinosWrappers { BlockSparseMatrix::BlockSparseMatrix () {} - - + + BlockSparseMatrix::~BlockSparseMatrix () {} @@ -94,13 +94,19 @@ namespace TrilinosWrappers const BlockSparsityPattern &block_sparsity_pattern) { const unsigned int n_block_rows = input_maps.size(); - - // Call the other reinit function ... + + Assert (n_block_rows == block_sparsity_pattern.n_block_rows(), + ExcDimensionMismatch (n_block_rows, + block_sparsity_pattern.n_block_rows())); + Assert (n_block_rows == block_sparsity_pattern.n_block_cols(), + ExcDimensionMismatch (n_block_rows, + block_sparsity_pattern.n_block_cols())); + + // Call the other basic reinit function ... reinit (n_block_rows, n_block_rows); - // ... and now assign the correct - // blocks. - + // ... and then assign the correct + // data to the blocks. for (unsigned int r=0; rn_block_rows(); ++r) for (unsigned int c=0; cn_block_cols(); ++c) { @@ -111,6 +117,37 @@ namespace TrilinosWrappers + void + BlockSparseMatrix:: + reinit (const std::vector &input_maps, + const ::dealii::BlockSparseMatrix &dealii_block_sparse_matrix, + const double drop_tolerance) + { + const unsigned int n_block_rows = input_maps.size(); + + Assert (n_block_rows == dealii_block_sparse_matrix.n_block_rows(), + ExcDimensionMismatch (n_block_rows, + dealii_block_sparse_matrix.n_block_rows())); + Assert (n_block_rows == dealii_block_sparse_matrix.n_block_cols(), + ExcDimensionMismatch (n_block_rows, + dealii_block_sparse_matrix.n_block_cols())); + + // Call the other basic reinit function ... + reinit (n_block_rows, n_block_rows); + + // ... and then assign the correct + // data to the blocks. + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + { + this->block(r,c).reinit(input_maps[r],input_maps[c], + dealii_block_sparse_matrix.block(r,c), + drop_tolerance); + } + } + + + void BlockSparseMatrix::compress() { diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index cd8838d66a..ef7c0ef741 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -26,23 +26,23 @@ namespace TrilinosWrappers SparseMatrix::const_iterator::Accessor:: visit_present_row () { - // if we are asked to visit the - // past-the-end line, then simply - // release all our caches and go on - // with life + // if we are asked to visit the + // past-the-end line, then simply + // release all our caches and go on + // with life if (this->a_row == matrix->m()) - { - colnum_cache.reset (); - value_cache.reset (); + { + colnum_cache.reset (); + value_cache.reset (); - return; - } + return; + } - // otherwise first flush Trilinos caches + // otherwise first flush Trilinos caches matrix->compress (); - // get a representation of the present - // row + // get a representation of the present + // row int ncols; int colnums = matrix->n(); TrilinosScalar *values = new TrilinosScalar(colnums); @@ -52,15 +52,15 @@ namespace TrilinosWrappers ncols, &(values[0])); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - // copy it into our caches if the line - // isn't empty. if it is, then we've - // done something wrong, since we - // shouldn't have initialized an - // iterator for an empty line (what - // would it point to?) + // copy it into our caches if the line + // isn't empty. if it is, then we've + // done something wrong, since we + // shouldn't have initialized an + // iterator for an empty line (what + // would it point to?) Assert (ncols != 0, ExcInternalError()); colnum_cache.reset (new std::vector (colnums, - colnums+ncols)); + colnums+ncols)); value_cache.reset (new std::vector (values, values+ncols)); } } @@ -68,17 +68,17 @@ namespace TrilinosWrappers // The constructor is actually the // only point where we have to check - // whether we build a serial or - // a parallel Trilinos matrix. - // In the end, it does not even - // matter how many threads there - // are, but only if we use an - // MPI compiler or a standard + // whether we build a serial or + // a parallel Trilinos matrix. + // In the end, it does not even + // matter how many threads there + // are, but only if we use an + // MPI compiler or a standard // compiler. So, one thread on - // an MPI compiler will still get - // a parallel interface. + // an MPI compiler will still get + // a parallel interface. SparseMatrix::SparseMatrix () - : + : #ifdef DEAL_II_COMPILER_SUPPORTS_MPI row_map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)), #else @@ -86,56 +86,56 @@ namespace TrilinosWrappers #endif col_map (row_map), matrix (std::auto_ptr - (new Epetra_FECrsMatrix(Copy, row_map, 0))), - last_action (Insert) + (new Epetra_FECrsMatrix(Copy, row_map, 0))), + last_action (Insert) {} SparseMatrix::SparseMatrix (const Epetra_Map &InputMap, const unsigned int n_max_entries_per_row) - : + : row_map (InputMap), col_map (row_map), matrix (std::auto_ptr - (new Epetra_FECrsMatrix(Copy, row_map, + (new Epetra_FECrsMatrix(Copy, row_map, int(n_max_entries_per_row), false))), - last_action (Insert) + last_action (Insert) {} SparseMatrix::SparseMatrix (const Epetra_Map &InputMap, const std::vector &n_entries_per_row) - : + : 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) + true))), + 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, + (new Epetra_FECrsMatrix(Copy, row_map, col_map, int(n_max_entries_per_row), false))), - last_action (Insert) + 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, row_map, col_map, (int*)const_cast(&(n_entries_per_row[0])), - false))), - last_action (Insert) + true))), + last_action (Insert) {} @@ -145,46 +145,6 @@ namespace TrilinosWrappers } - - void - SparseMatrix::reinit (const SparsityPattern &sparsity_pattern, - const unsigned int n_max_entries_per_row) - { - - 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; rowInsertGlobalValues(row, row_length, - &values[0], &row_indices[0]); - } - - // In the end, the matrix needs to - // be compressed in order to be - // really ready. However, that is - // a collective operation, so it - // has to be called on all processes - // by the user, whereas this function - // should only be used on one processor - // since our sparsity pattern data - // types are all serial. - } - - void SparseMatrix::reinit (const SparsityPattern &sparsity_pattern) @@ -200,6 +160,7 @@ namespace TrilinosWrappers // so do not check for consistent // column numbers here. // + // // this bug is filed in the Sandia // bugzilla under #4123 and should be // fixed for version 9.0 @@ -212,44 +173,40 @@ namespace TrilinosWrappers for (unsigned int row=0; row values; + std::vector row_indices; + + for (unsigned int row=0; rowInsertGlobalValues(row, row_length, + &values[0], &row_indices[0]); + } - reinit (sparsity_pattern, n_max_entries_per_row); + // In the end, the matrix needs to + // be compressed in order to be + // really ready. However, that is + // a collective operation, so it + // has to be called on all processes + // by the user, whereas this function + // should only be used on one processor + // since our sparsity pattern data + // types are all serial as of now. } void SparseMatrix::reinit (const Epetra_Map &input_map, - const SparsityPattern &sparsity_pattern) + const SparsityPattern &sparsity_pattern) { - matrix.reset(); - - unsigned int n_rows = sparsity_pattern.n_rows(); - - Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(), - ExcDimensionMismatch (input_map.NumGlobalElements(), - sparsity_pattern.n_rows())); - Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(), - ExcDimensionMismatch (input_map.NumGlobalElements(), - sparsity_pattern.n_cols())); - - row_map = input_map; - col_map = row_map; - - std::vector 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); + reinit (input_map, input_map, sparsity_pattern); } @@ -257,7 +214,7 @@ namespace TrilinosWrappers void SparseMatrix::reinit (const Epetra_Map &input_row_map, const Epetra_Map &input_col_map, - const SparsityPattern &sparsity_pattern) + const SparsityPattern &sparsity_pattern) { matrix.reset(); @@ -280,45 +237,53 @@ namespace TrilinosWrappers matrix = std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, col_map, - &n_entries_per_row[0], true)); + &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); + reinit (sparsity_pattern); } void SparseMatrix::reinit (const Epetra_Map &input_map, - const ::dealii::SparseMatrix &deal_ii_sparse_matrix, - const double drop_tolerance) + const ::dealii::SparseMatrix &dealii_sparse_matrix, + const double drop_tolerance) + { + reinit (input_map, input_map, dealii_sparse_matrix, drop_tolerance); + } + + + + void + SparseMatrix::reinit (const Epetra_Map &input_row_map, + const Epetra_Map &input_col_map, + const ::dealii::SparseMatrix &dealii_sparse_matrix, + const double drop_tolerance) { matrix.reset(); - unsigned int n_rows = deal_ii_sparse_matrix.m(); + unsigned int n_rows = dealii_sparse_matrix.m(); - Assert (input_map.NumGlobalElements() == (int)n_rows, - ExcDimensionMismatch (input_map.NumGlobalElements(), + Assert (input_row_map.NumGlobalElements() == (int)n_rows, + ExcDimensionMismatch (input_row_map.NumGlobalElements(), n_rows)); - Assert (input_map.NumGlobalElements() == (int)deal_ii_sparse_matrix.n(), - ExcDimensionMismatch (input_map.NumGlobalElements(), - deal_ii_sparse_matrix.n())); - - row_map = input_map; - col_map = row_map; + Assert (input_col_map.NumGlobalElements() == (int)dealii_sparse_matrix.n(), + ExcDimensionMismatch (input_col_map.NumGlobalElements(), + dealii_sparse_matrix.n())); + + row_map = input_row_map; + col_map = input_col_map; std::vector n_entries_per_row(n_rows); for (unsigned int row=0; row (new Epetra_FECrsMatrix(Copy, row_map, col_map, - &n_entries_per_row[0], true)); - + &n_entries_per_row[0], true)); + std::vector values; std::vector row_indices; @@ -329,8 +294,8 @@ namespace TrilinosWrappers unsigned int index = 0; for (dealii::SparseMatrix::const_iterator - p = deal_ii_sparse_matrix.begin(row); - p != deal_ii_sparse_matrix.end(row); ++p) + p = dealii_sparse_matrix.begin(row); + p != dealii_sparse_matrix.end(row); ++p) if (std::fabs(p->value()) > drop_tolerance) { row_indices[index] = p->column(); @@ -340,7 +305,7 @@ namespace TrilinosWrappers const int n_row_entries = index; matrix->InsertGlobalValues(row, n_row_entries, - &values[0], &row_indices[0]); + &values[0], &row_indices[0]); } } @@ -350,9 +315,9 @@ namespace TrilinosWrappers void SparseMatrix::clear () { - // When we clear the matrix, - // reset the pointer and - // generate an empty matrix. + // When we clear the matrix, + // reset the pointer and + // generate an empty matrix. matrix.reset(); #ifdef DEAL_II_COMPILER_SUPPORTS_MPI row_map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)), @@ -369,7 +334,7 @@ namespace TrilinosWrappers void SparseMatrix::compress () { - // flush buffers + // flush buffers int ierr; if (row_map.SameAs(col_map)) ierr = matrix->GlobalAssemble (); @@ -407,17 +372,17 @@ 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 == Add) { - int ierr; + int ierr; if (row_map.SameAs(col_map)) ierr = matrix->GlobalAssemble (false); else ierr = matrix->GlobalAssemble(col_map, row_map, false); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); last_action = Insert; } @@ -442,28 +407,28 @@ 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) { - int ierr; + int ierr; if (row_map.SameAs(col_map)) ierr = matrix->GlobalAssemble (false); else ierr = matrix->GlobalAssemble(col_map, row_map, false); - + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - last_action = Add; + last_action = Add; } - // we have to do above actions in any - // case to be consistent with the MPI - // communication model (see the - // comments in the documentation of - // TrilinosWrappers::Vector), but we - // can save some work if the addend is - // zero + // we have to do above actions in any + // case to be consistent with the MPI + // communication model (see the + // comments in the documentation of + // TrilinosWrappers::Vector), but we + // can save some work if the addend is + // zero if (value == 0) return; @@ -487,8 +452,8 @@ namespace TrilinosWrappers Assert (matrix->Filled()==true, ExcMessage("Matrix must be compressed before invoking clear_row.")); - // Only do this on the rows - // owned locally on this processor. + // Only do this on the rows + // owned locally on this processor. int local_row = matrix->LRID(row); if (local_row >= 0) { @@ -499,15 +464,15 @@ namespace TrilinosWrappers values, col_indices); Assert (ierr == 0, - ExcTrilinosError(ierr)); + ExcTrilinosError(ierr)); int* diag_find = std::find(col_indices,col_indices+num_entries, - local_row); + local_row); int diag_index = (int)(diag_find - col_indices); for (int j=0; j 0.) values[diag_index] = new_diag_value; @@ -530,17 +495,17 @@ namespace TrilinosWrappers SparseMatrix::el (const unsigned int i, const unsigned int j) const { - // Extract local indices in - // the matrix. + // Extract local indices in + // the matrix. int trilinos_i = matrix->LRID(i), trilinos_j = matrix->LRID(j); TrilinosScalar value = 0.; - // If the data is not on the - // present processor, we can't - // continue. + // If the data is not on the + // present processor, we can't + // 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 @@ -573,7 +538,7 @@ namespace TrilinosWrappers // finally get it. int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present, - trilinos_j); + trilinos_j); int el_index = (int)(el_find - col_indices); @@ -656,14 +621,14 @@ namespace TrilinosWrappers { Assert (row < m(), ExcInternalError()); - // get a representation of the present - // row + // get a representation of the present + // row int ncols = -1; int local_row = matrix->RowMap().LID(row); - // on the processor who owns this - // row, we'll have a non-negative - // value. + // on the processor who owns this + // row, we'll have a non-negative + // value. if (local_row >= 0) { int ierr = matrix->NumMyRowEntries (local_row, ncols); @@ -738,7 +703,7 @@ namespace TrilinosWrappers void SparseMatrix::vmult (Vector &dst, - const Vector &src) const + const Vector &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); @@ -768,7 +733,7 @@ namespace TrilinosWrappers void SparseMatrix::vmult_add (Vector &dst, - const Vector &src) const + const Vector &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); @@ -804,7 +769,7 @@ namespace TrilinosWrappers TrilinosScalar SparseMatrix::matrix_scalar_product (const Vector &u, - const Vector &v) const + const Vector &v) const { Vector tmp(v.map); vmult (tmp, v); @@ -828,10 +793,10 @@ namespace TrilinosWrappers // TODO: Currently this only flips - // a flag that tells Trilinos that - // any application should be done with - // the transpose. However, the matrix - // structure is not reset. + // a flag that tells Trilinos that + // any application should be done with + // the transpose. However, the matrix + // structure is not reset. void SparseMatrix::transpose () { @@ -839,13 +804,13 @@ namespace TrilinosWrappers if (!matrix->UseTranspose()) { - ierr = matrix->SetUseTranspose (true); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + ierr = matrix->SetUseTranspose (true); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } else { - ierr = matrix->SetUseTranspose (false); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + ierr = matrix->SetUseTranspose (false); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } } @@ -882,9 +847,9 @@ namespace TrilinosWrappers - // As of now, no particularly neat - // ouput is generated in case of - // multiple processors. + // As of now, no particularly neat + // ouput is generated in case of + // multiple processors. void SparseMatrix::print (std::ostream &out) const { double * values; @@ -893,7 +858,7 @@ namespace TrilinosWrappers for (int i=0; iNumMyRows(); ++i) { - matrix->ExtractMyRowView (i, num_entries, values, indices); + matrix->ExtractMyRowView (i, num_entries, values, indices); for (int j=0; jGRID(j)] << ") " << values[j] << std::endl; -- 2.39.5