From e0190b791100d9e1f95ce3043ed81d0024cc0fb1 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 10 Sep 2008 09:40:49 +0000 Subject: [PATCH] Some minor changes. git-svn-id: https://svn.dealii.org/trunk@16773 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/trilinos_sparse_matrix.h | 77 +++++++++------- deal.II/lac/source/trilinos_sparse_matrix.cc | 90 +++++++++++++++++-- 2 files changed, 130 insertions(+), 37 deletions(-) diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index ca2eaf08a8..1f037f175b 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -334,7 +334,9 @@ namespace TrilinosWrappers typedef TrilinosScalar value_type; /** - * Default constructor. + * Default + * constructor. Generates an + * empty (zero-size) matrix. */ SparseMatrix (); @@ -491,7 +493,7 @@ namespace TrilinosWrappers const Epetra_Map &input_col_map, const SparsityPattern &sparsity_pattern); - /** + /** * This function copies the * content in * sparse_matrix to @@ -642,11 +644,12 @@ namespace TrilinosWrappers * the matrix without having to * read entries (such as the * locations of non-zero - * elements) from it -- without - * this o peration, removing - * constraints on parallel - * matrices is a rather - * complicated procedure. + * elements) from it — + * without this operation, + * removing constraints on + * parallel matrices is a + * rather complicated + * procedure. * * The second parameter can be * used to set the diagonal @@ -683,20 +686,17 @@ namespace TrilinosWrappers * may be an expensive * operation and you should * always take care where to - * call this function. In - * contrast to the respective - * function in the @p - * SparseMatrix class, we don't - * throw an exception if the - * respective entry doesn't - * exist in the sparsity - * pattern of this class, since - * Trilinos does not transmit - * this information. On the - * other hand, an exception - * will be thrown when the - * requested element is not - * saved on the calling + * call this function. As in + * the &p SparseMatrix + * class, we throw an exception + * if the respective entry + * doesn't exist in the + * sparsity pattern of this + * class, which is requested + * from Trilinos. Moreover, an + * exception will be thrown + * when the requested element + * is not saved on the calling * process. * * This function is therefore @@ -1163,6 +1163,14 @@ namespace TrilinosWrappers << "An error with error number " << arg1 << " occured while calling a Trilinos function"); + /** + * Exception + */ + DeclException2 (ExcInvalidIndex, + int, int, + << "The entry with index <" << arg1 << ',' << arg2 + << "> does not exist."); + /** * Exception */ @@ -1248,6 +1256,9 @@ namespace TrilinosWrappers * elements. The actual type, * a sparse matrix, is set in * the constructor. + * + * TODO: This object should + * finally become private. */ std::auto_ptr matrix; @@ -1285,6 +1296,7 @@ namespace TrilinosWrappers } + inline unsigned int const_iterator::Accessor::column() const @@ -1294,6 +1306,7 @@ namespace TrilinosWrappers } + inline unsigned int const_iterator::Accessor::index() const @@ -1303,6 +1316,7 @@ namespace TrilinosWrappers } + inline TrilinosScalar const_iterator::Accessor::value() const @@ -1312,6 +1326,7 @@ namespace TrilinosWrappers } + inline const_iterator:: const_iterator(const SparseMatrix *matrix, @@ -1351,6 +1366,7 @@ namespace TrilinosWrappers } + inline const_iterator const_iterator::operator++ (int) @@ -1361,6 +1377,7 @@ namespace TrilinosWrappers } + inline const const_iterator::Accessor & const_iterator::operator* () const @@ -1369,6 +1386,7 @@ namespace TrilinosWrappers } + inline const const_iterator::Accessor * const_iterator::operator-> () const @@ -1377,6 +1395,7 @@ namespace TrilinosWrappers } + inline bool const_iterator:: @@ -1387,6 +1406,7 @@ namespace TrilinosWrappers } + inline bool const_iterator:: @@ -1396,6 +1416,7 @@ namespace TrilinosWrappers } + inline bool const_iterator:: @@ -1407,17 +1428,8 @@ namespace TrilinosWrappers } } - - - inline - TrilinosScalar - SparseMatrix::operator() (const unsigned int i, - const unsigned int j) const - { - return el(i,j); - } + - inline SparseMatrix::const_iterator @@ -1427,6 +1439,7 @@ namespace TrilinosWrappers } + inline SparseMatrix::const_iterator SparseMatrix::end() const @@ -1435,6 +1448,7 @@ namespace TrilinosWrappers } + inline SparseMatrix::const_iterator SparseMatrix::begin(const unsigned int r) const @@ -1447,6 +1461,7 @@ namespace TrilinosWrappers } + inline SparseMatrix::const_iterator SparseMatrix::end(const unsigned int r) const diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 0d00c4fb09..9488901c6b 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -376,9 +376,9 @@ namespace TrilinosWrappers // flush buffers int ierr; if (row_map.SameAs(col_map)) - ierr = matrix->GlobalAssemble (); + ierr = matrix->GlobalAssemble (true); else - ierr = matrix->GlobalAssemble (col_map, row_map); + ierr = matrix->GlobalAssemble (col_map, row_map, true); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); @@ -530,6 +530,82 @@ namespace TrilinosWrappers + TrilinosScalar + SparseMatrix::operator() (const unsigned int i, + const unsigned int j) const + { + // 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 throw + // an exception. This is on of + // the two tiny differences to + // the el(i,j) call, which does + // not throw any assertions. + if ((trilinos_i == -1 ) || (trilinos_j == -1)) + { + Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first, + local_range().second)); + } + else + { + // Check whether the matrix + // already is transformed to + // local indices. + if (matrix->Filled() == false) + matrix->FillComplete(true); + + // Prepare pointers for extraction + // of a view of the row. + int nnz_present = matrix->NumMyEntries(trilinos_i); + int nnz_extracted; + int *col_indices; + TrilinosScalar *values; + + // Generate the view and make + // sure that we have not generated + // an error. + int ierr = matrix->ExtractMyRowView(trilinos_i, nnz_extracted, + values, col_indices); + Assert (ierr==0, ExcTrilinosError(ierr)); + + Assert (nnz_present == nnz_extracted, + ExcDimensionMismatch(nnz_present, nnz_extracted)); + + // Search the index where we + // look for the value, and then + // finally get it. + + int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present, + trilinos_j); + + int el_index = (int)(el_find - col_indices); + + // This is actually the only + // difference to the el(i,j) + // function, which means that + // we throw an exception in + // this case instead of just + // returning zero for an + // element that is not present + // in the sparsity pattern. + if (!el_find) + { + Assert (false, ExcInvalidIndex (i,j)); + } + else + value = values[el_index]; + + } + + return value; + } + + + TrilinosScalar SparseMatrix::el (const unsigned int i, const unsigned int j) const @@ -544,7 +620,9 @@ namespace TrilinosWrappers // continue. Just print out // zero. - // TODO: Is this reasonable? + // TODO: Is this reasonable? Or + // should we retain the assert + // call? if ((trilinos_i == -1 ) || (trilinos_j == -1)) { return 0.; @@ -689,7 +767,7 @@ namespace TrilinosWrappers TrilinosScalar SparseMatrix::l1_norm () const { - if (!matrix->Filled()) + if (matrix->Filled() == false) matrix->FillComplete(); TrilinosScalar result = matrix->NormOne(); @@ -702,7 +780,7 @@ namespace TrilinosWrappers TrilinosScalar SparseMatrix::linfty_norm () const { - if (!matrix->Filled()) + if (matrix->Filled() == false) matrix->FillComplete(); TrilinosScalar result = matrix->NormInf(); @@ -715,7 +793,7 @@ namespace TrilinosWrappers TrilinosScalar SparseMatrix::frobenius_norm () const { - if (!matrix->Filled()) + if (matrix->Filled() == false) matrix->FillComplete(); TrilinosScalar result = matrix->NormFrobenius(); -- 2.39.5