From 0e0d9c3fe8bcf7fe36cabe7215b261f2e8f185b8 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 26 Aug 2008 15:36:33 +0000 Subject: [PATCH] Corrected some minor bugs in Trilinos wrappers. Added clear_row() functionality to sparse matrix class and a print() function in the deal.II style. git-svn-id: https://svn.dealii.org/trunk@16669 0785d39b-7218-0410-832d-ea1e28bc413d --- .../numerics/matrices.all_dimensions.cc | 17 +-- .../lac/include/lac/trilinos_sparse_matrix.h | 56 +++++--- deal.II/lac/include/lac/trilinos_vector.h | 9 +- .../lac/source/trilinos_precondition_amg.cc | 2 +- deal.II/lac/source/trilinos_sparse_matrix.cc | 105 ++++++++++++--- deal.II/lac/source/trilinos_vector.cc | 125 ++++++++++-------- 6 files changed, 204 insertions(+), 110 deletions(-) diff --git a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc index 78da534924..3034f81236 100644 --- a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc @@ -746,18 +746,11 @@ namespace TrilinosWrappers // then eliminate these rows and set // their diagonal entry to what we have - // determined above. note that for trilinos - // matrices interleaving read with write - // operations is very expensive. thus, we - // here always replace the diagonal - // element, rather than first checking - // whether it is nonzero and in that case - // preserving it. this is different from - // the case of deal.II sparse matrices - // treated in the other functions. -//TODO: clear_row is not currently implemented for Trilinos - Assert (false, ExcInternalError()); -// matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry); + // determined above. if the value already + // is nonzero, it will be preserved, + // in accordance with the basic + // matrix classes in deal.II. + matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry); // the next thing is to set right hand // side to the wanted value. there's one diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 3285cc6a63..8d3bdfd5d4 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -477,6 +477,21 @@ namespace TrilinosWrappers */ void clear (); + /** + * Trilinos matrices store their own + * sparsity patterns. So, in analogy to + * our own SparsityPattern class, + * this function compresses the + * sparsity pattern and allows the + * resulting matrix to be used in all + * other operations where before only + * assembly functions were + * allowed. This function must + * therefore be called once you have + * assembled the matrix. + */ + void compress (); + /** * This operator assigns a scalar to a * matrix. Since this does usually not @@ -545,7 +560,10 @@ namespace TrilinosWrappers * pattern, though (but retains * the allocated memory in case * new entries are again added - * later). + * later). Note that this + * is a global operation, so this + * needs to be done on all + * MPI processes. * * This operation is used in * eliminating constraints (e.g. due to @@ -554,7 +572,7 @@ namespace TrilinosWrappers * the matrix without having to read * entries (such as the locations of * non-zero elements) from it -- - * without this operation, removing + * without this o peration, removing * constraints on parallel matrices is * a rather complicated procedure. * @@ -563,8 +581,8 @@ namespace TrilinosWrappers * to a value different from zero. The * default is to set it to zero. */ - // void clear_row (const unsigned int row, - // const TrilinosScalar new_diag_value = 0); + void clear_row (const unsigned int row, + const TrilinosScalar new_diag_value = 0); /** * Same as clear_row(), except that it @@ -579,24 +597,9 @@ namespace TrilinosWrappers * the diagonal entries, you have to * set them by hand. */ - // void clear_rows (const std::vector &rows, - // const TrilinosScalar new_diag_value = 0); - - /** - * Trilinos matrices store their own - * sparsity patterns. So, in analogy to - * our own SparsityPattern class, - * this function compresses the - * sparsity pattern and allows the - * resulting matrix to be used in all - * other operations where before only - * assembly functions were - * allowed. This function must - * therefore be called once you have - * assembled the matrix. - */ - void compress (); - + void clear_rows (const std::vector &rows, + const TrilinosScalar new_diag_value = 0); + /** * Return the value of the entry * (i,j). This may be an @@ -998,6 +1001,15 @@ namespace TrilinosWrappers */ void write_ascii (); + /** + * Print the matrix to the given + * stream, using the format + * (line,col) value, + * i.e. one nonzero entry of the + * matrix per line. + */ + void print (std::ostream &out) const; + // TODO: Write an overloading // of the operator << for output. // Since the underlying Trilinos diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h index 7b9e151a6d..33af1a0c38 100755 --- a/deal.II/lac/include/lac/trilinos_vector.h +++ b/deal.II/lac/include/lac/trilinos_vector.h @@ -770,6 +770,13 @@ namespace TrilinosWrappers void ratio (const Vector &a, const Vector &b); + /** + * Output of vector in user-defined + * format in analogy to the + * dealii::Vector class. + */ + void print (const char* format = 0) const; + /** * Print to a * stream. @p precision denotes @@ -900,7 +907,7 @@ namespace TrilinosWrappers * * @ingroup TrilinosWrappers * @relates TrilinosWrappers::Vector - * @author Wolfgang Bangerth, 2004 + * @author Martin Kronbichler, Wolfgang Bangerth, 2008 */ inline void swap (Vector &u, Vector &v) diff --git a/deal.II/lac/source/trilinos_precondition_amg.cc b/deal.II/lac/source/trilinos_precondition_amg.cc index 5791c8a2fc..7d287c230d 100755 --- a/deal.II/lac/source/trilinos_precondition_amg.cc +++ b/deal.II/lac/source/trilinos_precondition_amg.cc @@ -175,7 +175,7 @@ namespace TrilinosWrappers Epetra_Vector LHS (View, multigrid_operator->OperatorDomainMap(), dst.begin()); - Epetra_Vector RHS (View, multigrid_operator->OperatorDomainMap(), + Epetra_Vector RHS (View, multigrid_operator->OperatorRangeMap(), const_cast(src.begin())); const int res = multigrid_operator->ApplyInverse (RHS, LHS); diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index d0b2054fd5..3af6a7a2cf 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -195,14 +195,10 @@ namespace TrilinosWrappers ExcDimensionMismatch (matrix->NumGlobalRows(), sparsity_pattern.n_rows())); - // Trilinos seems to have a bug for - // rectangular matrices at this point, - // 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 + // Trilinos seems to have a bug for + // rectangular matrices at this point, + // so do not check for consistent + // column numbers here. // Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(), // ExcDimensionMismatch (matrix->NumGlobalCols(), // sparsity_pattern.n_cols())); @@ -331,7 +327,7 @@ namespace TrilinosWrappers for (dealii::SparseMatrix::const_iterator p = deal_ii_sparse_matrix.begin(row); p != deal_ii_sparse_matrix.end(row); ++p) - if (std::abs(p->value()) > drop_tolerance) + if (std::fabs(p->value()) > drop_tolerance) { row_indices[index] = p->column(); values[index] = p->value(); @@ -461,7 +457,7 @@ namespace TrilinosWrappers // case to be consistent with the MPI // communication model (see the // comments in the documentation of - // TrilinosWrappers::MPI::Vector), but we + // TrilinosWrappers::Vector), but we // can save some work if the addend is // zero if (value == 0) @@ -480,6 +476,52 @@ namespace TrilinosWrappers + void + SparseMatrix::clear_row (const unsigned int row, + const TrilinosScalar new_diag_value) + { + Assert (matrix->Filled()==true, + ExcMessage("Matrix must be compressed before invoking clear_row.")); + + // Only do this on the rows + // owned locally on this processor. + int local_row = matrix->LRID(row); + if (local_row >= 0) + { + TrilinosScalar *values; + int *col_indices; + int num_entries; + const int ierr = matrix->ExtractMyRowView(local_row, num_entries, + values, col_indices); + + Assert (ierr == 0, + ExcTrilinosError(ierr)); + + int* diag_find = std::find(col_indices,col_indices+num_entries, + local_row); + int diag_index = (int)(diag_find - col_indices); + + for (int j=0; j 0.) + values[diag_index] = new_diag_value; + } + } + + + + void + SparseMatrix::clear_rows (const std::vector &rows, + const TrilinosScalar new_diag_value) + { + for (unsigned int row=0; rowLRID(i), trilinos_j = matrix->LRID(j); TrilinosScalar value = 0.; - + // If the data is not on the // present processor, we can't // continue. @@ -525,16 +567,17 @@ namespace TrilinosWrappers // Search the index where we // look for the value, and then // finally get it. - int* index = std::find(&col_indices[0],&col_indices[0] + nnz_present, - trilinos_j); + + int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present, + trilinos_j); + + int el_index = (int)(el_find - col_indices); - int position; - if (!index) + if (!el_find) value = 0; else { - position = (int)(index - &(col_indices[0])); - value = values[position]; + value = values[el_index]; } } @@ -548,8 +591,6 @@ namespace TrilinosWrappers { Assert (m() == n(), ExcNotQuadratic()); - // this doesn't seem to work any - // different than any other element return el(i,i); } @@ -589,7 +630,7 @@ namespace TrilinosWrappers { int begin, end; begin = matrix->RowMap().MinMyGID(); - end = matrix->RowMap().MaxMyGID(); + end = matrix->RowMap().MaxMyGID()+1; return std::make_pair (begin, end); } @@ -827,12 +868,36 @@ namespace TrilinosWrappers return false; } + + void SparseMatrix::write_ascii () { Assert (false, ExcNotImplemented()); } + + + // As of now, no particularly neat + // ouput is generated in case of + // multiple processors. + void SparseMatrix::print (std::ostream &out) const + { + double * values; + int * indices; + int num_entries; + + for (int i=0; iNumMyRows(); ++i) + { + matrix->ExtractMyRowView (i, num_entries, values, indices); + for (int j=0; jGRID(j)] << ") " + << values[j] << std::endl; + } + + AssertThrow (out, ExcIO()); + } + } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc index 0c952be353..09d6bfd169 100755 --- a/deal.II/lac/source/trilinos_vector.cc +++ b/deal.II/lac/source/trilinos_vector.cc @@ -207,7 +207,7 @@ namespace TrilinosWrappers { int begin, end; begin = vector->Map().MinMyGID(); - end = vector->Map().MaxMyGID(); + end = vector->Map().MaxMyGID()+1; return std::make_pair (begin, end); } @@ -520,17 +520,14 @@ namespace TrilinosWrappers ExcMessage("The given value is not finite but " "either infinite or Not A Number (NaN)")); - std::vector list (size(), s); - - int* index = new int[size()]; - for (unsigned int i=0; iSumIntoGlobalValues(size(), index, &list[0]); + unsigned int n_local = local_size(); + int ierr; - delete[] index; - - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + for (unsigned int i=0; iSumIntoMyValue(i,0,s); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } } @@ -545,14 +542,14 @@ namespace TrilinosWrappers void Vector::add (const TrilinosScalar a, - const Vector &v) + const Vector &v) { Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but " "either infinite or Not A Number (NaN)")); - const int ierr = vector->Update(a, *(v.vector), 1); + const int ierr = vector->Update(a, *(v.vector), 1.); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -560,9 +557,9 @@ namespace TrilinosWrappers void Vector::add (const TrilinosScalar a, - const Vector &v, - const TrilinosScalar b, - const Vector &w) + const Vector &v, + const TrilinosScalar b, + const Vector &w) { Assert (numbers::is_finite(a), @@ -572,7 +569,7 @@ namespace TrilinosWrappers ExcMessage("The given value is not finite but " "either infinite or Not A Number (NaN)")); - const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.0); + const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -581,14 +578,14 @@ namespace TrilinosWrappers void Vector::sadd (const TrilinosScalar s, - const Vector &v) + const Vector &v) { Assert (numbers::is_finite(s), ExcMessage("The given value is not finite but " "either infinite or Not A Number (NaN)")); - const int ierr = vector->Update(1.0, *(v.vector), s); + const int ierr = vector->Update(1., *(v.vector), s); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -597,8 +594,8 @@ namespace TrilinosWrappers void Vector::sadd (const TrilinosScalar s, - const TrilinosScalar a, - const Vector &v) + const TrilinosScalar a, + const Vector &v) { Assert (numbers::is_finite(s), @@ -617,10 +614,10 @@ namespace TrilinosWrappers void Vector::sadd (const TrilinosScalar s, - const TrilinosScalar a, - const Vector &v, - const TrilinosScalar b, - const Vector &w) + const TrilinosScalar a, + const Vector &v, + const TrilinosScalar b, + const Vector &w) { Assert (numbers::is_finite(s), @@ -642,12 +639,12 @@ namespace TrilinosWrappers void Vector::sadd (const TrilinosScalar s, - const TrilinosScalar a, - const Vector &v, - const TrilinosScalar b, - const Vector &w, - const TrilinosScalar c, - const Vector &x) + const TrilinosScalar a, + const Vector &v, + const TrilinosScalar b, + const Vector &w, + const TrilinosScalar c, + const Vector &x) { Assert (numbers::is_finite(s), @@ -663,12 +660,13 @@ namespace TrilinosWrappers ExcMessage("The given value is not finite but " "either infinite or Not A Number (NaN)")); - // Update member can only input two other vectors so + // Update member can only input + // two other vectors so // do it in two steps const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - - const int jerr = vector->Update(c, *(x.vector), 1.0); + + const int jerr = vector->Update(c, *(x.vector), 1.); AssertThrow (jerr == 0, ExcTrilinosError(jerr)); } @@ -692,11 +690,10 @@ namespace TrilinosWrappers ExcMessage("The given value is not finite but " "either infinite or Not A Number (NaN)")); - Assert (size() == v.size(), - ExcDimensionMismatch (size(), v.size())); - - const int ierr = vector->Update(a, *(v.vector), 0.0); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + *vector = *v.vector; + map = v.map; + + *this *= a; } @@ -715,15 +712,12 @@ namespace TrilinosWrappers ExcMessage("The given value is not finite but " "either infinite or Not A Number (NaN)")); - Assert (size() == v.size(), - ExcDimensionMismatch (size(), v.size())); - - Assert (size() == w.size(), - ExcDimensionMismatch (size(), w.size())); - - const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 0.0); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + Assert (v.size() == w.size(), + ExcDimensionMismatch (v.size(), w.size())); + *vector = *v.vector; + map = v.map; + sadd (a, b, w); } @@ -732,8 +726,8 @@ namespace TrilinosWrappers Vector::ratio (const Vector &v, const Vector &w) { - Assert (size() == v.size(), - ExcDimensionMismatch (size(), v.size())); + Assert (v.size() == w.size(), + ExcDimensionMismatch (v.size(), w.size())); Assert (size() == w.size(), ExcDimensionMismatch (size(), w.size())); @@ -745,6 +739,28 @@ namespace TrilinosWrappers } + + // TODO: up to now only local data + // printed out! Find a way to neatly + // output distributed data... + void + Vector::print (const char *format) const + { + Assert (vector->GlobalLength()!=0, ExcEmptyObject()); + + for (unsigned int j=0; jExtractView (&val, &leading_dimension); @@ -790,10 +807,10 @@ namespace TrilinosWrappers // Just swap the pointers to the // two Epetra vectors that hold all // the data. - std::auto_ptr tmp; - tmp = v.vector; - v.vector = vector; - vector = tmp; + Vector *p_v = &v, *p_this = this; + Vector* tmp = p_v; + p_v = p_this; + p_this = tmp; } -- 2.39.5