From: kronbichler Date: Thu, 6 Nov 2008 09:52:15 +0000 (+0000) Subject: The Trilinos sparse matrix class can now insert elements using set(i,j) as long as... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=997c1afa55f9221b83aa5772b116b29e653a03a1;p=dealii-svn.git The Trilinos sparse matrix class can now insert elements using set(i,j) as long as compress() is not called. git-svn-id: https://svn.dealii.org/trunk@17483 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 6520591907..b557e10e76 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -836,15 +836,17 @@ namespace TrilinosWrappers * Set the element (i,j) * to @p value. * - * Just as the respective call in - * deal.II SparseMatrix - * class (but in contrast to the - * situation for PETSc based - * matrices), this function - * throws an exception if an - * entry does not exist in the - * sparsity pattern. Moreover, if - * value is not a finite + * This function is able to insert + * new elements into the matrix as + * long as compress() has not been + * called, so the sparsity pattern + * will be extended. When compress() + * is called for the first time, then + * this is no longer possible and an + * insertion of elements at positions + * which have not been initialized + * will throw an exception. Moreover, + * if value is not a finite * number an exception is thrown. */ void set (const unsigned int i, diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 8feeac441b..dfc1be2f26 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -630,9 +630,21 @@ namespace TrilinosWrappers int trilinos_i = i; int trilinos_j = j; - const int ierr = matrix->ReplaceGlobalValues (trilinos_i, 1, - const_cast(&value), - &trilinos_j); + int ierr; + + // If the matrix is not yet filled, + // we can insert new entries into + // the matrix using this + // 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); + else + ierr = matrix->ReplaceGlobalValues (trilinos_i, 1, + const_cast(&value), + &trilinos_j); AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j)); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); @@ -762,7 +774,7 @@ namespace TrilinosWrappers // already is transformed to // local indices. if (matrix->Filled() == false) - matrix->FillComplete(col_map, row_map, true); + matrix->GlobalAssemble(col_map, row_map, true); // Prepare pointers for extraction // of a view of the row. @@ -785,10 +797,10 @@ namespace TrilinosWrappers // look for the value, and then // finally get it. - int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present, + int* el_find = std::find(col_indices, col_indices + nnz_present, trilinos_j); - - int el_index = (int)(el_find - col_indices); + + int local_col_index = (int)(el_find - col_indices); // This is actually the only // difference to the el(i,j) @@ -798,13 +810,12 @@ namespace TrilinosWrappers // returning zero for an // element that is not present // in the sparsity pattern. - if (!el_find) + if (local_col_index == nnz_present) { Assert (false, ExcInvalidIndex (i,j)); } else - value = values[el_index]; - + value = values[local_col_index]; } return value; @@ -841,7 +852,7 @@ namespace TrilinosWrappers // already is transformed to // local indices. if (!matrix->Filled()) - matrix->FillComplete(col_map, row_map, true); + matrix->GlobalAssemble(col_map, row_map, true); // Prepare pointers for extraction // of a view of the row. @@ -863,18 +874,23 @@ namespace TrilinosWrappers // 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); + int* el_find = std::find(col_indices, col_indices + nnz_present, + trilinos_j); + + int local_col_index = (int)(el_find - col_indices); - if (!el_find) + + // This is actually the only + // difference to the () function + // querying (i,j), where we throw an + // exception instead of just + // returning zero for an element + // that is not present in the + // sparsity pattern. + if (local_col_index == nnz_present) value = 0; else - { - value = values[el_index]; - } + value = values[local_col_index]; } return value;