From: Martin Kronbichler Date: Fri, 27 Aug 2010 14:26:03 +0000 (+0000) Subject: Cleanup reinit functions. Remove some unnecessary code, make the code scale. X-Git-Tag: v8.0.0~5648 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5984b4e8e22e05d35e50eba6b389b07a89519987;p=dealii.git Cleanup reinit functions. Remove some unnecessary code, make the code scale. git-svn-id: https://svn.dealii.org/trunk@21750 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 7efa621178..7e6814ae9d 100644 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -373,25 +373,43 @@ namespace TrilinosWrappers temp_vector.clear(); matrix.reset(); + // if we want to exchange data, build + // a usual Trilinos sparsity pattern + // and let that handle the + // exchange. otherwise, manually + // create a CrsGraph, which consumes + // considerably less memory because it + // can set correct number of indices + // right from the start + if (exchange_data) + { + SparsityPattern trilinos_sparsity; + trilinos_sparsity.reinit (input_row_map, input_col_map, + sparsity_pattern, exchange_data); + reinit (trilinos_sparsity); + + return; + } + const unsigned int n_rows = sparsity_pattern.n_rows(); Assert (exchange_data == false, ExcNotImplemented()); if (input_row_map.Comm().MyPID() == 0) { - Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(), - ExcDimensionMismatch (input_row_map.NumGlobalElements(), - sparsity_pattern.n_rows())); - Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(), - ExcDimensionMismatch (input_col_map.NumGlobalElements(), - sparsity_pattern.n_cols())); + AssertDimension (sparsity_pattern.n_rows(), + static_cast(input_row_map.NumGlobalElements())); + AssertDimension (sparsity_pattern.n_cols(), + static_cast(input_col_map.NumGlobalElements())); } column_space_map.reset (new Epetra_Map (input_col_map)); - std::vector n_entries_per_row(n_rows); + const unsigned int first_row = input_row_map.MinMyGID(), + last_row = input_row_map.MaxMyGID()+1; + std::vector n_entries_per_row(last_entry - first_entry); - for (unsigned int row=0; row graph; - int * first_row_length = &n_entries_per_row[input_row_map.MinMyGID()]; if (input_row_map.Comm().NumProc() > 1) - graph.reset (new Epetra_CrsGraph (Copy, input_row_map, first_row_length, - true)); + graph.reset (new Epetra_CrsGraph (Copy, input_row_map, + &n_entries_per_row[0], true)); else graph.reset (new Epetra_CrsGraph (Copy, input_row_map, input_col_map, - first_row_length, true)); + &n_entries_per_row[0], true)); // This functions assumes that the // sparsity pattern sits on all processors // (completely). The parallel version uses // an Epetra graph that is already // distributed. - Assert (graph->NumGlobalRows() == (int)sparsity_pattern.n_rows(), - ExcDimensionMismatch (graph->NumGlobalRows(), - sparsity_pattern.n_rows())); - // now insert the indices + // now insert the indices std::vector row_indices; - for (unsigned int row=0; rowEpetra_CrsGraph::InsertGlobalIndices (row, row_length, - &row_indices[0]); - } + typename SparsityType::row_iterator col_num = sp.row_begin (row), + row_end = sp.row_end(row); + for (unsigned int col = 0; col_num != row_end; ++col_num, ++col) + row_indices[col] = *col_num; + + graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length, + &row_indices[0]); + } // Eventually, optimize the graph // structure (sort indices, make memory @@ -461,13 +476,11 @@ namespace TrilinosWrappers // check whether we got the number of // columns right. - Assert (graph->NumGlobalCols() == (int)sparsity_pattern.n_cols(), - ExcDimensionMismatch (graph->NumGlobalCols(), - sparsity_pattern.n_cols())); + AssertDimension (sparsity_pattern.n_cols(), + static_cast(graph->NumGlobalCols())); // And now finally generate the matrix. matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false)); - last_action = Zero; // In the end, the matrix needs to @@ -478,30 +491,16 @@ namespace TrilinosWrappers - template <> - void - SparseMatrix::reinit (const Epetra_Map &input_row_map, - const Epetra_Map &input_col_map, - const CompressedSimpleSparsityPattern &sparsity_pattern, - const bool exchange_data) - { - SparsityPattern trilinos_sparsity; - trilinos_sparsity.reinit (input_row_map, input_col_map, sparsity_pattern, - exchange_data); - reinit (trilinos_sparsity); - } - - - void SparseMatrix::reinit (const SparsityPattern &sparsity_pattern) { + temp_vector.clear (); + matrix.reset (); + // reinit with a (parallel) Trilinos // sparsity pattern. column_space_map.reset (new Epetra_Map (sparsity_pattern.domain_partitioner())); - temp_vector.clear (); - matrix.reset (); matrix.reset (new Epetra_FECrsMatrix (Copy, sparsity_pattern.trilinos_sparsity_pattern(), false)); compress(); diff --git a/deal.II/lac/source/trilinos_sparsity_pattern.cc b/deal.II/lac/source/trilinos_sparsity_pattern.cc index 9c7a70ec29..43775e5aad 100644 --- a/deal.II/lac/source/trilinos_sparsity_pattern.cc +++ b/deal.II/lac/source/trilinos_sparsity_pattern.cc @@ -242,10 +242,8 @@ namespace TrilinosWrappers { // release memory before reallocation graph.reset (); - Assert (n_entries_per_row.size() == - static_cast(input_row_map.NumGlobalElements()), - ExcDimensionMismatch (n_entries_per_row.size(), - input_row_map.NumGlobalElements())); + AssertDimension (n_entries_per_row.size(), + static_cast(input_row_map.NumGlobalElements())); column_space_map.reset (new Epetra_Map (input_col_map)); compressed = false; @@ -282,14 +280,10 @@ namespace TrilinosWrappers { graph.reset (); - Assert (sp.n_rows() == - static_cast(input_row_map.NumGlobalElements()), - ExcDimensionMismatch (sp.n_rows(), - input_row_map.NumGlobalElements())); - Assert (sp.n_cols() == - static_cast(input_col_map.NumGlobalElements()), - ExcDimensionMismatch (sp.n_cols(), - input_col_map.NumGlobalElements())); + AssertDimension (sp.n_rows(), + static_cast(input_row_map.NumGlobalElements())); + AssertDimension (sp.n_cols(), + static_cast(input_col_map.NumGlobalElements())); column_space_map.reset (new Epetra_Map (input_col_map)); compressed = false; @@ -297,35 +291,32 @@ namespace TrilinosWrappers Assert (input_row_map.LinearMap() == true, ExcMessage ("This function is not efficient if the map is not contiguous.")); - std::vector n_entries_per_row(input_row_map.MaxMyGID()- - input_row_map.MinMyGID() + 1); + const unsigned int first_row = input_row_map.MinMyGID(), + last_row = input_row_map.MaxMyGID()+1; + std::vector n_entries_per_row(last_entry - first_entry); - for (unsigned int row=input_row_map.MinMyGID(); - row(input_row_map.MaxMyGID()+1); - ++row) - n_entries_per_row[row-input_row_map.MinMyGID()] = sp.row_length(row); + for (unsigned int row=first_row; row 1) graph.reset(new Epetra_FECrsGraph(Copy, input_row_map, - n_entries_per_row[0], + &n_entries_per_row[0], false)); else graph.reset (new Epetra_FECrsGraph(Copy, input_row_map, input_col_map, - n_entries_per_row[0], + &n_entries_per_row[0], false)); - Assert (graph->NumGlobalRows() == (int)sp.n_rows(), - ExcDimensionMismatch (graph->NumGlobalRows(), - sp.n_rows())); + AssertDimension (sp.n_rows(), + static_cast(graph->NumGlobalRows())); - const unsigned int n_rows = sp.n_rows(); std::vector row_indices; // Include possibility to exchange data // since CompressedSimpleSparsityPattern is // able to do so - for (unsigned int row=0; rowEpetra_CrsGraph::InsertGlobalIndices (row, row_length, &row_indices[0]); } - else if ( exchange_data ) + else + for (unsigned int row=0; rowInsertGlobalIndices (1, (int*)&row, row_length, &row_indices[0]); + graph->InsertGlobalIndices (1, static_cast(&row), row_length, + &row_indices[0]); } compress();