From: kronbichler Date: Fri, 21 Nov 2008 10:01:43 +0000 (+0000) Subject: Added the way elements are inserted/replaced/added into Trilinos matrices. Adding... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=84f1e68cae997c4a785d50548e5f0309300fc962;p=dealii-svn.git Added the way elements are inserted/replaced/added into Trilinos matrices. Adding elements into the global stiffness matrix one by one using add() should now take between 30 and 50 precent less time than it did before. git-svn-id: https://svn.dealii.org/trunk@17674 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 36f23e6009..7f63ba7ade 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -1243,8 +1243,8 @@ void BoussinesqFlowProblem::solve () GridTools::minimal_cell_diameter(triangulation) / std::max (get_maximal_velocity(), 1.e-5); - std::cout << " " << "Time step: " << time_step - << std::endl; + pcout << " " << "Time step: " << time_step + << std::endl; temperature_solution = old_temperature_solution; diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index b04db9f9e8..e06d4ef268 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -734,17 +734,13 @@ namespace TrilinosWrappers void compress (); /** - * Returns the state of the - * matrix, i.e., whether - * compress() needs to be called - * after an operation requiring - * data exchange. Does only - * return non-true values when - * used in debug mode, - * since it is quite expensive to - * keep track of all operations - * that lead to the need for - * compress(). + * Returns the state of the matrix, + * i.e., whether compress() needs to + * be called after an operation + * requiring data exchange. A call to + * compress() is also after the + * method set() is called (even when + * working in serial). */ bool is_compressed () const; //@} diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 054c80a60e..6b33306a8f 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -208,6 +208,8 @@ namespace TrilinosWrappers SparseMatrix & SparseMatrix::copy_from (const SparseMatrix &m) { + row_map = m.row_map; + col_map = m.col_map; *matrix = *m.matrix; return *this; } @@ -313,8 +315,8 @@ namespace TrilinosWrappers // ExcDimensionMismatch (matrix->NumGlobalCols(), // sparsity_pattern.n_cols())); - std::vector values; - std::vector row_indices; + std::vector values; + std::vector row_indices; for (unsigned int row=0; rowInsertGlobalValues (row, row_length, - &values[0], &row_indices[0]); + matrix->Epetra_CrsMatrix::InsertGlobalValues ((int)row, row_length, + &values[0], + &row_indices[0]); } last_action = Zero; @@ -412,7 +415,7 @@ namespace TrilinosWrappers // ExcDimensionMismatch (matrix->NumGlobalCols(), // sparsity_pattern.n_cols())); - std::vector values; + std::vector values; std::vector row_indices; for (unsigned int row=0; rowInsertGlobalValues (row, row_length, - &values[0], &row_indices[0]); + matrix->Epetra_CrsMatrix::InsertGlobalValues ((int)row, row_length, + &values[0], + &row_indices[0]); } last_action = Zero; @@ -531,7 +535,7 @@ namespace TrilinosWrappers (new Epetra_FECrsMatrix(Copy, row_map, &n_entries_per_row[0], false)); - std::vector values; + std::vector values; std::vector row_indices; for (unsigned int row=0; rowvalue(); ++index; } - const int n_row_entries = index; - matrix->InsertGlobalValues(row, n_row_entries, - &values[0], &row_indices[0]); + const int row_length = index; + + matrix->Epetra_CrsMatrix::InsertGlobalValues ((int)row, row_length, + &values[0], + &row_indices[0]); } compress(); @@ -680,38 +686,80 @@ namespace TrilinosWrappers if (last_action != Insert) last_action = Insert; -#ifdef DEBUG - if (in_local_range (i) == false) - compressed = false; -#endif - - int trilinos_i = i; - int trilinos_j = j; + compressed = false; int ierr; + // In case this matrix owns the + // actual row, we can save some + // computing time by directly + // accessing the Epetra_CrsMatrix + // set functions. + if (row_map.MyGID(i) == true) + { // 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); - if (ierr > 0) - ierr = 0; + if (matrix->Filled() == false) + { + ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues ( + (int)i, 1, + const_cast(&value), + (int*)&j); + + // When adding up elements, we do + // not want to create exceptions in + // the case when adding elements. + if (ierr > 0) + ierr = 0; + } + else + { + ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues ( + (int)i, 1, + const_cast(&value), + (int*)&j); + } } else { - ierr = matrix->ReplaceGlobalValues (trilinos_i, 1, - const_cast(&value), - &trilinos_j); + int trilinos_i = i; + int trilinos_j = j; + const TrilinosScalar* value_ptr = &value; + + + // 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 (1, &trilinos_i, + 1, &trilinos_j, + &value_ptr, + Epetra_FECrsMatrix::ROW_MAJOR); + + // When adding up elements, we do + // not want to create exceptions in + // the case when adding elements. + if (ierr > 0) + ierr = 0; + } + else + { + ierr = matrix->ReplaceGlobalValues (1, &trilinos_i, + 1, &trilinos_j, + &value_ptr, + Epetra_FECrsMatrix::ROW_MAJOR); + } } - AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j)); + Assert (ierr <= 0, ExcAccessToNonPresentElement(i,j)); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } @@ -747,17 +795,43 @@ namespace TrilinosWrappers if (value == 0) return; -#ifdef DEBUG - if (in_local_range (i) == false) - compressed = false; -#endif + int ierr; + // If the calling matrix owns the row to + // which we want to add this value, we + // can directly call the Epetra_CrsMatrix + // input function, which is much faster + // than the Epetra_FECrsMatrix function. + if (row_map.MyGID(i) == true) + { + ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues( + (int)i, 1, + const_cast(&value), + (int*)&j); + } + else + { + // When we're at off-processor data, we + // have to stick with the standard + // SumIntoGlobalValues + // function. Nevertheless, the way we + // call it is the fastest one (any other + // will lead to repeated allocation and + // deallocation of memory in order to + // call the function we already use, + // which is very unefficient if writing + // one element at a time). + compressed = false; - int trilinos_i = i; - int trilinos_j = j; + int trilinos_i = i; + int trilinos_j = j; - const int ierr = matrix->SumIntoGlobalValues (trilinos_i, 1, - const_cast(&value), - &trilinos_j); + const TrilinosScalar* value_ptr = &value; + + ierr = matrix->SumIntoGlobalValues (1, &trilinos_i, + 1, &trilinos_j, + &value_ptr, + Epetra_FECrsMatrix::ROW_MAJOR); + } AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j)); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); @@ -1262,26 +1336,53 @@ namespace TrilinosWrappers = std::max (max_row_length, static_cast(rhs.matrix->NumGlobalEntries(row))); - std::vector column_indices (max_row_length); - std::vector values (max_row_length); - - for (unsigned int row=local_range.first; - row < local_range.second; ++row) + std::vector column_indices (max_row_length); + std::vector values (max_row_length); + + // If both matrices have been transformed + // to local index space (in Trilinos + // speak: they are filled), we can + // extract MyRows instead of GlobalRows, + // which is faster. + if (matrix->Filled() == true && rhs.matrix->Filled() == true) + for (unsigned int row=local_range.first; + row < local_range.second; ++row) + { + const int row_local = matrix->RowMap().LID(row); + int n_entries; + rhs.matrix->ExtractMyRowCopy (row_local, max_row_length, + n_entries, + &values[0], + &column_indices[0]); + for (int i=0; iSumIntoMyValues (row_local, n_entries, value_ptr, + &column_indices[0]); + } + else { - int n_entries; - rhs.matrix->ExtractGlobalRowCopy (row, max_row_length, - n_entries, - &values[0], - &column_indices[0]); - for (int i=0; iSumIntoGlobalValues (row, n_entries, - &values[0], - &column_indices[0]); + for (unsigned int row=local_range.first; + row < local_range.second; ++row) + { + int n_entries; + rhs.matrix->Epetra_CrsMatrix::ExtractGlobalRowCopy ((int)row, + max_row_length, + n_entries, + &values[0], + &column_indices[0]); + for (int i=0; iEpetra_CrsMatrix::SumIntoGlobalValues ((int)row, + n_entries, + &values[0], + &column_indices[0]); + } + compress (); } - - compress (); }